
!  +++++++++++++++++++++++ BOXTETRA_ROUTINES ++++++++++++++++++++++++++

MODULE BOXTETRA_ROUTINES
USE PRECISION_PARAMETERS
USE MATH_FUNCTIONS, ONLY: CROSS_PRODUCT
IMPLICIT NONE

PRIVATE

INTEGER, DIMENSION(0:3,0:5), TARGET :: BOX_PLANE2VERT
INTEGER, DIMENSION(0:3,0:5), TARGET :: BOX_PLANE2EDGE
INTEGER, DIMENSION(0:1,0:11), TARGET :: BOX_EDGE2VERT

INTEGER, DIMENSION(0:2,0:3), TARGET :: TETRA_PLANE2VERT
INTEGER, DIMENSION(0:2,0:3), TARGET :: TETRA_PLANE2EDGE
INTEGER, DIMENSION(0:1,0:5) :: TETRA_EDGE2VERT

REAL(EB), DIMENSION(0:2,0:5), TARGET :: BOX_NORMALS=RESHAPE((/&
  -1.0_EB,  0.0_EB,  0.0_EB,&
   1.0_EB,  0.0_EB,  0.0_EB,&
   0.0_EB,-1.0_EB,   0.0_EB,&
   0.0_EB,  1.0_EB,  0.0_EB,&
   0.0_EB,  0.0_EB, -1.0_EB,&
   0.0_EB,  0.0_EB,  1.0_EB&
   /),(/3,6/))

REAL(EB), DIMENSION(0:2,0:3), TARGET :: TETRA_NORMALS, TETRA_VERTS
REAL(EB), DIMENSION(0:2,0:7), TARGET :: BOX_VERTS

REAL(EB), DIMENSION(0:5) :: TETRA_BOUNDS
REAL(EB), DIMENSION(0:5,0:3),TARGET :: TETRA_PLANE_BOUNDS
REAL(EB), DIMENSION(0:5,0:5),TARGET :: BOX_PLANE_BOUNDS

INTEGER, PARAMETER :: MIN_X=0, MAX_X=1, MIN_Y=2, MAX_Y=3, MIN_Z=4, MAX_Z=5

INTEGER :: I, J

!       6-----------7
!      /.          /|
!    /  .        /  |
!   4-----------5   |
!   |   .       |   |
!   |   .       |   |
!   |   2.......|...3
!   |  .        |  /
!   | .         | /
!   |.          |/
!   0-----------1

! BOX_PLANE2VERT(edge,plane)
DATA ( (BOX_PLANE2VERT(I,J), I=0,3),J=0,5) /&
  0,2,4,6,&
  1,3,5,7,&
  0,1,4,5,&
  2,3,6,7,&
  0,1,2,3,&
  4,5,6,7 &
  /

!       ------7------
!      /.           /
!     2 .         3 |
!    /  .        /  |
!   ------6------   |
!   |  10       |  11
!   |   .       |   |
!   |   .....5..|...|
!   8  .        9  /
!   | 0         | 1
!   |.          |/
!   |----4------|

DATA ( (BOX_PLANE2EDGE(I,J), I=0,3),J=0,5) /&
  0,2,8,10,&
  1,3,9,11,&
  4,6,8,9,&
  5,7,10,11,&
  0,1,4,5,&
  2,3,6,7 &
  /

!       6-----7-----7
!      /.           /
!     2 .         3 |
!    /  .        /  |
!   4-----6-----5   |
!   |  10       |  11
!   |   .       |   |
!   |   2....5..|...3
!   8  .        9  /
!   | 0         | 1
!   |.          |/
!   0----4------1
! planes: 0-left 1-right 2-front 3-back 4-bottom 5-top
! edges: 0-bottom left  1-bottom right 2-top left   3-top right
!        4-bottom front 5-bottom back  6-top front  7-top back
!        8-front left   9-front right 10-back left 11-back right
! vertices: 0-bottom left front 1-bottom right front 2-bottom left back 3-bottom right back
!           4-top left front    5-top right front    6-top left back    7-top right back

 DATA ( (BOX_EDGE2VERT(I,J), I=0,1), J=0,11) /&
  0,2,  1,3,  4,6,  5,7,&
  0,1,  2,3,  4,5,  6,7,&
  0,4,  1,5,  2,6,  3,7&
  /

!           3
!          /.\
!         / . \
!        /  5  \
!       /   .   \
!      3    2    4
!     /   .   .   \
!    /  2       1  \
!   / .           . \
!  0-------0---------1

DATA ( (TETRA_PLANE2VERT(I,J), I=0,2),J=0,3) /&
  0,3,1,&
  1,3,2,&
  0,2,3,&
  0,1,2&
  /

DATA ( (TETRA_PLANE2EDGE(I,J), I=0,2),J=0,3) /&
  0,3,4,&
  1,4,5,&
  2,5,3,& !double check (was 2 3 5)
  0,1,2&
  /

 DATA ( (TETRA_EDGE2VERT(I,J), I=0,1), J=0,5) /&
  0,1,  1,2,  2,0,&
  0,3,  1,3,  2,3&
  /

PUBLIC GET_TETRABOX_VOLUME, GET_VERTS, TETRAHEDRON_VOLUME, REMOVE_DUPLICATE_VERTS, GET_POLYHEDRON_VOLUME, TEST_IN_TETRA0,&
       GEOMCLIP, DECIMATE, DECIMATE_FB, IN_TRIANGLE, IN_POLYGON, IS_ANGLE_GE_180, POLY2TRI

CONTAINS

!  ------------------ GET_TETRABOX_VOLUME ------------------------

SUBROUTINE GET_TETRABOX_VOLUME(BOX_BOUNDS,V0,V1,V2,V3,TETRABOX_VOLUME,AREAS,CENTROID)

! compute the volume of the intersection of a box and a tetrahedron

REAL(EB), DIMENSION(0:5), INTENT(IN) :: BOX_BOUNDS
REAL(EB), DIMENSION(0:2), INTENT(IN) :: V0, V1, V2, V3
REAL(EB), INTENT(OUT) :: TETRABOX_VOLUME,AREAS(6), CENTROID(3)

REAL(EB), DIMENSION(0:599) :: VERTS
INTEGER :: NVERTS, NSIDES, NFACES
INTEGER, DIMENSION(0:599) :: FACES
INTEGER, DIMENSION(0:199) :: WHICH_SIDE
INTEGER :: SIDE_INDEX(0:9)
INTEGER :: BOX_STATE(-1:5)
INTEGER :: VERT_STATE

AREAS=0.0_EB
CENTROID=0.0_EB

! obtain vertices of intersection region

CALL GET_VERTS(BOX_BOUNDS,V0,V1,V2,V3,VERTS,NVERTS,FACES,SIDE_INDEX,WHICH_SIDE,NFACES,NSIDES,&
               BOX_STATE,VERT_STATE,AREAS)

IF (VERT_STATE==1) THEN  ! box is inside tetrahedron
   TETRABOX_VOLUME = (BOX_BOUNDS(1)-BOX_BOUNDS(0))*(BOX_BOUNDS(3)-BOX_BOUNDS(2))*(BOX_BOUNDS(5)-BOX_BOUNDS(4))
   CENTROID(1)=(BOX_BOUNDS(0)+BOX_BOUNDS(1))/2.0_EB
   CENTROID(2)=(BOX_BOUNDS(2)+BOX_BOUNDS(3))/2.0_EB
   CENTROID(3)=(BOX_BOUNDS(4)+BOX_BOUNDS(5))/2.0_EB
   RETURN
ELSE IF (VERT_STATE==2) THEN ! tetrahedron is inside box
   TETRABOX_VOLUME = ABS(TETRAHEDRON_VOLUME(V0,V1,V2,V3))
   CENTROID(1:3) = (V0(0:2)+V1(0:2)+V2(0:2)+V3(0:2))/4.0_EB
   RETURN
ELSE IF (VERT_STATE==3) THEN ! tetrahedron and box do not intersect
   TETRABOX_VOLUME = 0.0_EB
   CENTROID=0.0_EB
   RETURN
ENDIF

! compute volume using computed vertices

CALL GET_POLYHEDRON_VOLUME(VERTS,NVERTS,FACES,NFACES,TETRABOX_VOLUME,CENTROID)

RETURN
END SUBROUTINE GET_TETRABOX_VOLUME

!  ------------------ GET_AREAS ------------------------

SUBROUTINE GET_AREAS(VERTS,FACES,NFACES,WHICH_SIDE,AREAS)
INTEGER, INTENT(IN) :: NFACES
REAL(EB), DIMENSION(0:599), INTENT(IN), TARGET :: VERTS
INTEGER, DIMENSION(0:599),TARGET, INTENT(IN) :: FACES
INTEGER, DIMENSION(0:199), INTENT(IN) :: WHICH_SIDE
REAL(EB), INTENT(OUT) :: AREAS(6)

INTEGER :: I
REAL(EB), POINTER, DIMENSION(:) :: V1,V2,V3
INTEGER, POINTER,DIMENSION(:) :: F
REAL(EB) :: AREA
INTEGER :: SIDE

AREAS=0.0_EB
DO I = 0, NFACES-1
   F(0:2)=>FACES(3*I:3*I+2)
   V1(1:3)=>VERTS(3*F(0):3*F(0)+2)
   V2(1:3)=>VERTS(3*F(1):3*F(1)+2)
   V3(1:3)=>VERTS(3*F(2):3*F(2)+2)
   AREA = TRIANGLE_AREA(V1,V2,V3)
   SIDE = WHICH_SIDE(I) + 1
   IF (SIDE<=6)AREAS(SIDE) = AREAS(SIDE) + AREA
ENDDO

END SUBROUTINE GET_AREAS

!  ------------------ TETRAHEDRON_VOLUME ------------------------

!              D1
!             /|\
!            / | \
!           /  |  \
!          /   |   \
!         /    |    \
!        /     B4    \
!       /     . .     \
!      /     .    .    \
!     /    .        .   \
!    /   .            .  \
!   /  .               .  \
!  / .                    .\
! C2------------------------A3

REAL(EB) FUNCTION TETRAHEDRON_VOLUME(A,B,C,D)

! determine the volume of a tetrahedron formed from vertices A, B, C and D

REAL(EB), DIMENSION(0:2), INTENT(IN) :: A, B, C, D
REAL(EB), DIMENSION(0:2) :: AMC, BMC, DMC, ACROSSB

AMC = A - C
BMC = B - C
DMC = D - C
CALL CROSS_PRODUCT(ACROSSB,AMC, BMC)
TETRAHEDRON_VOLUME = DOT_PRODUCT(ACROSSB,DMC)/6.0_EB
END FUNCTION TETRAHEDRON_VOLUME

!  ------------------ TRIANGLE_AREA ------------------------

!              C
!             / \
!            /   \
!           /     \
!          /       \
!         /         \
!        /           \
!       /             \
!      /               \
!     /                 \
!    /                   \
!   /                     \
!  /                       \
! A-------------------------B

REAL(EB) FUNCTION TRIANGLE_AREA(A,B,C)

! determine the area of a triangle formed from vertices A, B, and C

REAL(EB), INTENT(IN) :: A(3), B(3), C(3)
REAL(EB) :: AMC(3), BMC(3), ACROSSB(3)

AMC = A - C
BMC = B - C
CALL CROSS_PRODUCT(ACROSSB,AMC, BMC)
TRIANGLE_AREA = SQRT(DOT_PRODUCT(ACROSSB,ACROSSB))/2.0_EB
END FUNCTION TRIANGLE_AREA

!  ------------------ GET_POLYHEDRON_VOLUME ------------------------

SUBROUTINE GET_POLYHEDRON_VOLUME(VERTS,NVERTS,FACES,NFACES,POLYHEDRON_VOLUME,CENTROID)

REAL(EB), DIMENSION(0:599), INTENT(IN), TARGET :: VERTS
INTEGER, INTENT(IN) :: NVERTS, NFACES
INTEGER, TARGET, DIMENSION(0:599), INTENT(IN) :: FACES
REAL(EB), INTENT(OUT) :: POLYHEDRON_VOLUME, CENTROID(3)

REAL(EB), DIMENSION(0:2) :: V_CENTER
REAL(EB), DIMENSION(:), POINTER :: V0, V1, V2
INTEGER, POINTER, DIMENSION(:) :: FACEPTR
REAL(EB) :: VOL_TET

! compute center of intersection region

CENTROID(1:3)=0.0_EB
POLYHEDRON_VOLUME = 0.0_EB
IF (NVERTS==0) RETURN

V_CENTER(0:2)=0.0_EB
DO I = 0, NVERTS - 1
   V_CENTER(0:2) = V_CENTER(0:2) + VERTS(3*I:3*I+2)
ENDDO
V_CENTER = V_CENTER/REAL(NVERTS,EB)

! sum volumes of each region formed using polyhedron face and center

DO I = 0, NFACES-1
   FACEPTR(0:2) => FACES(3*I:3*I+2)
   V0(0:2) => VERTS(3*FACEPTR(0):3*FACEPTR(0)+2)
   V1(0:2) => VERTS(3*FACEPTR(1):3*FACEPTR(1)+2)
   V2(0:2) => VERTS(3*FACEPTR(2):3*FACEPTR(2)+2)
   VOL_TET = ABS(TETRAHEDRON_VOLUME(V0,V1,V2,V_CENTER))
   CENTROID(1:3) = CENTROID(1:3) + VOL_TET*(V0(0:2)+V1(0:2)+V2(0:2)+V_CENTER(0:2))/4.0_EB
   POLYHEDRON_VOLUME = POLYHEDRON_VOLUME + VOL_TET
ENDDO

IF (POLYHEDRON_VOLUME>0.0_EB) THEN
   CENTROID(1:3) = CENTROID(1:3)/POLYHEDRON_VOLUME
ELSE
   CENTROID(1:3) = 0.0_EB
ENDIF

END SUBROUTINE GET_POLYHEDRON_VOLUME

!  ------------------ GET_VERTS ------------------------

SUBROUTINE GET_VERTS(BOX_BOUNDS,VV0,VV1,VV2,VV3,VERTS,NVERTS,FACES,SIDE_INDEX,WHICH_SIDE,NFACES,NSIDES,&
                     BOX_STATE,VERT_STATE,AREAS)

! determine vertices of tetrahedron box intersection region ordered by faces

REAL(EB), DIMENSION(0:5), INTENT(IN) :: BOX_BOUNDS
REAL(EB), DIMENSION(0:2), INTENT(IN) :: VV0, VV1, VV2, VV3
REAL(EB), DIMENSION(0:599), INTENT(OUT), TARGET :: VERTS
INTEGER, DIMENSION(0:599),TARGET, INTENT(OUT) :: FACES
INTEGER, DIMENSION(0:199), INTENT(OUT) :: WHICH_SIDE
INTEGER, DIMENSION(0:9), INTENT(OUT) :: SIDE_INDEX
REAL(EB), INTENT(OUT) :: AREAS(6)

INTEGER, INTENT(OUT) :: NVERTS,NSIDES,NFACES,VERT_STATE
INTEGER, INTENT(OUT) :: BOX_STATE(-1:5)

INTEGER, DIMENSION(0:199) :: SIDE_START, SIDE_N

REAL(EB), DIMENSION(0:2) :: VERT
REAL(EB), POINTER, DIMENSION(:) :: BOXVERT, BOXNORMAL,  BOXVERT0, BOXVERT1
REAL(EB), POINTER, DIMENSION(:) :: TETRAVERT, TETRANORMAL, TETRAVERT0, TETRAVERT1
REAL(EB), POINTER, DIMENSION(:) :: FACEVERTS
INTEGER, POINTER, DIMENSION(:) :: EDGE
REAL(EB), DIMENSION(:), POINTER :: B_PLANE_BOUNDS, T_PLANE_BOUNDS
INTEGER :: TETRA_STATE(0:3)
REAL(EB) :: VERTS2D(0:299)
INTEGER :: NVERTS2D
INTEGER :: VOLUS(1)
INTEGER, PARAMETER :: MAXVERTS=600,MAXFACES=600,MAXVOLUS=1,NVOLUS=0
INTEGER :: FIRST_FACE_INDEX
REAL(EB), POINTER, DIMENSION(:) :: VPTR
INTEGER, POINTER, DIMENSION(:) :: FPTR
INTEGER, PARAMETER :: DEBUG=0
REAL(EB), PARAMETER :: EPS_SINGLE=1.0E-4_EB

INTEGER :: V,E,F,BP,BE,BV,TP,TE,TV

REAL(EB), DIMENSION(0:2) :: V0, V1, V2, V3

REAL(EB), POINTER, DIMENSION(:) :: TET0, TET1, TET2

! setup data structures

V0=VV0
V1=VV1
V2=VV2
V3=VV3
IF (TETRAHEDRON_VOLUME(V0,V1,V2,V3)<0.0_EB) THEN
! swap vertices if orientation is reversed
   V0=VV0
   V1=VV2
   V2=VV1
   V3=VV3
ENDIF
CALL SETUP_VERTS(BOX_BOUNDS,V0,V1,V2,V3)

NVERTS=0
NSIDES=0
SIDE_START = 0
SIDE_N = 0
SIDE_INDEX = 0
AREAS(1:6)=0.0_EB

! return if there is not an intersection

VERT_STATE = 3
IF (BOX_BOUNDS(1)<TETRA_BOUNDS(0) .OR. BOX_BOUNDS(0)>TETRA_BOUNDS(1)) RETURN
IF (BOX_BOUNDS(3)<TETRA_BOUNDS(2) .OR. BOX_BOUNDS(2)>TETRA_BOUNDS(3)) RETURN
IF (BOX_BOUNDS(5)<TETRA_BOUNDS(4) .OR. BOX_BOUNDS(4)>TETRA_BOUNDS(5)) RETURN

VERT_STATE=1
DO I = 0, 7
   IF (IN_TETRA(BOX_VERTS(0:2,I),-1)==0) THEN
      VERT_STATE=0
      EXIT
   ENDIF
ENDDO
IF (VERT_STATE==1) THEN ! box is inside tetrahedron
! define verts and faces for a box
   VERTS(0:2)=BOX_VERTS(0:2,0)
   VERTS(3:5)=BOX_VERTS(0:2,1)
   VERTS(6:8)=BOX_VERTS(0:2,2)
   VERTS(9:11)=BOX_VERTS(0:2,3)
   VERTS(12:14)=BOX_VERTS(0:2,4)
   VERTS(15:17)=BOX_VERTS(0:2,5)
   VERTS(18:20)=BOX_VERTS(0:2,6)
   VERTS(21:23)=BOX_VERTS(0:2,7)
   NVERTS=8
   FACES(0:5)  =(/0,1,5,0,5,4/)
   FACES(6:11) =(/1,3,7,1,7,5/)
   FACES(12:17)=(/3,2,6,3,6,7/)
   FACES(18:23)=(/2,0,4,2,4,6/)
   FACES(24:29)=(/4,5,7,4,7,6/)
   FACES(30:35)=(/2,3,1,2,1,0/)
   NFACES=12
   WHICH_SIDE(0:11) = (/2,2,1,1,3,3,0,0,5,5,4,4/)

! find area of each box side that intersects with the tetrahedron

   CALL GET_AREAS(VERTS,FACES,NFACES,WHICH_SIDE,AREAS)
   RETURN
ENDIF

IF (IN_BOX(V0,-1)==1 .AND. IN_BOX(V1,-1)==1 .AND. IN_BOX(V2,-1)==1 .AND. IN_BOX(V3,-1)==1) THEN ! tetrahedron is inside box
! define verts and faces for a tetrahedron
   VERTS(0:2)=TETRA_VERTS(0:2,0)
   VERTS(3:5)=TETRA_VERTS(0:2,1)
   VERTS(6:8)=TETRA_VERTS(0:2,2)
   VERTS(9:11)=TETRA_VERTS(0:2,3)
   NVERTS=4
   FACES(0:2)=(/0,1,3/)
   FACES(3:5)=(/1,2,3/)
   FACES(6:8)=(/0,3,2/)
   FACES(9:11)=(/0,2,1/)
   NFACES=4
   WHICH_SIDE(0:3) = (/6,7,8,9/)
   VERT_STATE=2

! find area of each box side that intersects with the tetrahedron

   CALL GET_AREAS(VERTS,FACES,NFACES,WHICH_SIDE,AREAS)
   RETURN
ENDIF

VERT_STATE=0

! determine vertices of intersection region

BOX_STATE = -1

! TETRA_STATE(I)=J (with J>=0) means that tetra plane I is in the same plane as box plane J
TETRA_STATE(0) = ON_BOX_PLANE(BOX_BOUNDS,V0,V3,V1)
TETRA_STATE(1) = ON_BOX_PLANE(BOX_BOUNDS,V1,V3,V2)
TETRA_STATE(2) = ON_BOX_PLANE(BOX_BOUNDS,V0,V2,V3)
TETRA_STATE(3) = ON_BOX_PLANE(BOX_BOUNDS,V0,V1,V2)

! BOX_STATE(I)=J (with J>=0) means that box plane I is in the same plane as tetra plane J
BOX_STATE(TETRA_STATE(0))=0
BOX_STATE(TETRA_STATE(1))=1
BOX_STATE(TETRA_STATE(2))=2
BOX_STATE(TETRA_STATE(3))=3

! for each box plane ...

DO BP = 0, 5
   V = BOX_PLANE2VERT(0,BP)
   BOXVERT(0:2) => BOX_VERTS(0:2,V)
   BOXNORMAL(0:2) => BOX_NORMALS(0:2,BP)
   B_PLANE_BOUNDS(0:5) => BOX_PLANE_BOUNDS(0:5,BP)

   IF (BOX_STATE(BP)/=-1 ) THEN ! this box plane coincides with a tetrahedon plane
                                  ! use a 2D method to find intersections
      TET0(0:2) => TETRA_VERTS(0:2,TETRA_PLANE2VERT(0,BOX_STATE(BP)))
      TET1(0:2) => TETRA_VERTS(0:2,TETRA_PLANE2VERT(1,BOX_STATE(BP)))
      TET2(0:2) => TETRA_VERTS(0:2,TETRA_PLANE2VERT(2,BOX_STATE(BP)))
      CALL GET_VERTS2D(B_PLANE_BOUNDS,BP,TET0,TET2,TET1,VERTS2D,NVERTS2D)
      IF (NVERTS2D>0) THEN
         VERTS(3*NVERTS:3*NVERTS+3*NVERTS2D-1)=VERTS2D(0:3*NVERTS2D-1)
         NVERTS = NVERTS + NVERTS2D
         SIDE_N(NSIDES) = SIDE_N(NSIDES) + NVERTS2D
      ENDIF
      IF (SIDE_N(NSIDES)>0) THEN
         SIDE_INDEX(NSIDES) = BP
         NSIDES = NSIDES + 1
      ENDIF
      CYCLE
   ENDIF


!  add intersection of tetrahedron edge and box plane (if on box)

   DO TE = 0, 5 ! tetra edge
      TETRAVERT0(0:2) => TETRA_VERTS(0:2,TETRA_EDGE2VERT(0,TE))
      TETRAVERT1(0:2) => TETRA_VERTS(0:2,TETRA_EDGE2VERT(1,TE))
      IF (B_PLANE_BOUNDS(0)>MAX(TETRAVERT0(0),TETRAVERT1(0))+EPS_SINGLE) CYCLE
      IF (B_PLANE_BOUNDS(1)<MIN(TETRAVERT0(0),TETRAVERT1(0))-EPS_SINGLE) CYCLE
      IF (B_PLANE_BOUNDS(2)>MAX(TETRAVERT0(1),TETRAVERT1(1))+EPS_SINGLE) CYCLE
      IF (B_PLANE_BOUNDS(3)<MIN(TETRAVERT0(1),TETRAVERT1(1))-EPS_SINGLE) CYCLE
      IF (B_PLANE_BOUNDS(4)>MAX(TETRAVERT0(2),TETRAVERT1(2))+EPS_SINGLE) CYCLE
      IF (B_PLANE_BOUNDS(5)<MIN(TETRAVERT0(2),TETRAVERT1(2))-EPS_SINGLE) CYCLE
      IF ( BOXPLANE_EDGE_INTERSECTION( BOXVERT,BP,TETRAVERT0,TETRAVERT1,VERT)/=1) CYCLE
      IF ( IN_BOX(VERT, BP)/=1) CYCLE

      VERTS(3*NVERTS:3*NVERTS+2)=VERT(0:2)
      SIDE_N(NSIDES) = SIDE_N(NSIDES) + 1
      NVERTS=NVERTS+1
   ENDDO

!  add intersection of box plane edge and tetrahedron plane (if on tetrahedron)

   EDGE(0:3) => BOX_PLANE2EDGE(0:3,BP)
   DO BE = 0, 3 ! box plane edges
      E = EDGE(BE)
      BOXVERT0(0:2) => BOX_VERTS(0:2,BOX_EDGE2VERT(0,E))
      BOXVERT1(0:2) => BOX_VERTS(0:2,BOX_EDGE2VERT(1,E))
      DO TP = 0, 3 ! tetra plane
         T_PLANE_BOUNDS(0:5) => TETRA_PLANE_BOUNDS(0:5,TP)
         IF (T_PLANE_BOUNDS(0)>MAX(BOXVERT0(0),BOXVERT1(0))+EPS_SINGLE) CYCLE
         IF (T_PLANE_BOUNDS(1)<MIN(BOXVERT0(0),BOXVERT1(0))-EPS_SINGLE) CYCLE
         IF (T_PLANE_BOUNDS(2)>MAX(BOXVERT0(1),BOXVERT1(1))+EPS_SINGLE) CYCLE
         IF (T_PLANE_BOUNDS(3)<MIN(BOXVERT0(1),BOXVERT1(1))-EPS_SINGLE) CYCLE
         IF (T_PLANE_BOUNDS(4)>MAX(BOXVERT0(2),BOXVERT1(2))+EPS_SINGLE) CYCLE
         IF (T_PLANE_BOUNDS(5)<MIN(BOXVERT0(2),BOXVERT1(2))-EPS_SINGLE) CYCLE
         V = TETRA_PLANE2VERT(0,TP)
         TETRAVERT(0:2) => TETRA_VERTS(0:2,V)
         TETRANORMAL(0:2) => TETRA_NORMALS(0:2,TP)
         IF ( PLANE_EDGE_INTERSECTION( TETRAVERT,TETRANORMAL,BOXVERT0,BOXVERT1,VERT)==1) THEN
            IF (IN_TETRA(VERT, TP)==1) THEN
               VERTS(3*NVERTS:3*NVERTS+2)=VERT(0:2)
               SIDE_N(NSIDES) = SIDE_N(NSIDES) + 1
               NVERTS=NVERTS+1
            ENDIF
         ENDIF
      ENDDO
   ENDDO

! add box plane vertices if inside tetrahedron

   DO BV = 0, 3
      BOXVERT(0:2) => BOX_VERTS(0:2,BOX_PLANE2VERT(BV,BP))
      IF (IN_TETRA(BOXVERT, -1)==1) THEN
         VERTS(3*NVERTS:3*NVERTS+2)=BOXVERT(0:2)
         SIDE_N(NSIDES) = SIDE_N(NSIDES) + 1
         NVERTS=NVERTS+1
      ENDIF
   ENDDO

   IF (SIDE_N(NSIDES)>0) THEN
      SIDE_INDEX(NSIDES) = BP
      NSIDES = NSIDES + 1
   ENDIF

ENDDO

! for each tetrahedron plane ....

DO TP = 0, 3
   IF (TETRA_STATE(TP)/=-1 ) CYCLE ! this tetrahedron plane coincides with a box plane
                                   ! in this case, use a 2D method to find intersections

   V = TETRA_PLANE2VERT(0,TP)
   TETRAVERT(0:2) => TETRA_VERTS(0:2,V)
   TETRANORMAL(0:2) => TETRA_NORMALS(0:2,TP)
   T_PLANE_BOUNDS(0:5) => TETRA_PLANE_BOUNDS(0:5,TP)

!  add intersection of box edge and tetrahedron plane (if on tetrahedron)

   DO BE = 0, 11 ! box edge
      BOXVERT0(0:2) => BOX_VERTS(0:2,BOX_EDGE2VERT(0,BE))
      BOXVERT1(0:2) => BOX_VERTS(0:2,BOX_EDGE2VERT(1,BE))
      IF (T_PLANE_BOUNDS(0)>MAX(BOXVERT0(0),BOXVERT1(0))+EPS_SINGLE) CYCLE
      IF (T_PLANE_BOUNDS(1)<MIN(BOXVERT0(0),BOXVERT1(0))-EPS_SINGLE) CYCLE
      IF (T_PLANE_BOUNDS(2)>MAX(BOXVERT0(1),BOXVERT1(1))+EPS_SINGLE) CYCLE
      IF (T_PLANE_BOUNDS(3)<MIN(BOXVERT0(1),BOXVERT1(1))-EPS_SINGLE) CYCLE
      IF (T_PLANE_BOUNDS(4)>MAX(BOXVERT0(2),BOXVERT1(2))+EPS_SINGLE) CYCLE
      IF (T_PLANE_BOUNDS(5)<MIN(BOXVERT0(2),BOXVERT1(2))-EPS_SINGLE) CYCLE
      IF ( PLANE_EDGE_INTERSECTION( TETRAVERT,TETRANORMAL,BOXVERT0,BOXVERT1,VERT)/=1) CYCLE
      IF (IN_TETRA(VERT, TP)/=1) CYCLE

      VERTS(3*NVERTS:3*NVERTS+2)=VERT(0:2)
      SIDE_N(NSIDES) = SIDE_N(NSIDES) + 1
      NVERTS=NVERTS+1
   ENDDO

!  add intersection of tetrahedron plane edge and box plane (if on box)

   EDGE(0:2) => TETRA_PLANE2EDGE(0:2,TP)
   DO TE = 0, 2 ! tetrahedron plane edges
      E = EDGE(TE)
      TETRAVERT0(0:2) => TETRA_VERTS(0:2,TETRA_EDGE2VERT(0,E))
      TETRAVERT1(0:2) => TETRA_VERTS(0:2,TETRA_EDGE2VERT(1,E))
      DO BP = 0, 5 ! box plane
         B_PLANE_BOUNDS(0:5) => BOX_PLANE_BOUNDS(0:5,BP)
         IF (B_PLANE_BOUNDS(0)>MAX(TETRAVERT0(0),TETRAVERT1(0))+EPS_SINGLE) CYCLE
         IF (B_PLANE_BOUNDS(1)<MIN(TETRAVERT0(0),TETRAVERT1(0))-EPS_SINGLE) CYCLE
         IF (B_PLANE_BOUNDS(2)>MAX(TETRAVERT0(1),TETRAVERT1(1))+EPS_SINGLE) CYCLE
         IF (B_PLANE_BOUNDS(3)<MIN(TETRAVERT0(1),TETRAVERT1(1))-EPS_SINGLE) CYCLE
         IF (B_PLANE_BOUNDS(4)>MAX(TETRAVERT0(2),TETRAVERT1(2))+EPS_SINGLE) CYCLE
         IF (B_PLANE_BOUNDS(5)<MIN(TETRAVERT0(2),TETRAVERT1(2))-EPS_SINGLE) CYCLE
         V = BOX_PLANE2VERT(0,BP)
         BOXVERT(0:2) => BOX_VERTS(0:2,V)
         BOXNORMAL(0:2) => BOX_NORMALS(0:2,BP)
         IF ( BOXPLANE_EDGE_INTERSECTION( BOXVERT,BP,TETRAVERT0,TETRAVERT1,VERT)/=1) CYCLE
         IF (IN_BOX(VERT, BP)/=1) CYCLE

         VERTS(3*NVERTS:3*NVERTS+2)=VERT(0:2)
         SIDE_N(NSIDES) = SIDE_N(NSIDES) + 1
         NVERTS=NVERTS+1
      ENDDO
   ENDDO

! add tetrahedron plane vertices if inside box

   DO TV = 0, 2
      TETRAVERT(0:2) => TETRA_VERTS(0:2,TETRA_PLANE2VERT(TV,TP))
      IF (IN_BOX(TETRAVERT, -1)==1) THEN
         VERTS(3*NVERTS:3*NVERTS+2)=TETRAVERT(0:2)
         SIDE_N(NSIDES) = SIDE_N(NSIDES) + 1
         NVERTS=NVERTS+1
      ENDIF
   ENDDO

   IF (SIDE_N(NSIDES)>0) THEN
      SIDE_INDEX(NSIDES) = 6+TP
      NSIDES = NSIDES + 1
   ENDIF

ENDDO

! determine vertex index at the start of each face

DO F = 1, NSIDES
  SIDE_START(F) = SIDE_START(F-1) + SIDE_N(F-1)
ENDDO

! order vertices of each face clockwise

DO F = 0, NSIDES-1
  IF (SIDE_N(F)<=3) CYCLE
  FACEVERTS(0:3*SIDE_N(F)-1) => VERTS(3*SIDE_START(F):3*SIDE_START(F+1)-1)
  CALL ORDER_VERTS(FACEVERTS,SIDE_N(F),SIDE_INDEX(F))
ENDDO

NFACES=0
DO I = 0, NSIDES-1
   IF (SIDE_N(I)<3) CYCLE
   DO J = SIDE_START(I)+1, SIDE_START(I+1)-2
      FACES(3*NFACES:3*NFACES+2) = (/SIDE_START(I),J,J+1/)
      WHICH_SIDE(NFACES) = SIDE_INDEX(I)
      NFACES = NFACES + 1
   ENDDO
ENDDO

FIRST_FACE_INDEX=0
CALL REMOVE_DUPLICATE_VERTS(NVERTS,NFACES,NVOLUS,MAXVERTS,MAXFACES,MAXVOLUS,FIRST_FACE_INDEX, &
                            VERTS,FACES,VOLUS,EPS_SINGLE)

CALL REMOVE_DEGENERATE_FACES(NFACES,MAXFACES,FACES,WHICH_SIDE)
IF (DEBUG==1) THEN
  DO I = 0, NVERTS-1
     VPTR(0:2) => VERTS(3*I:3*I+2)
     WRITE(6,*)I,VPTR(0),VPTR(1),VPTR(2)
  ENDDO
  DO I = 0, NFACES-1
     FPTR(0:2) => FACES(3*I:3*I+2)
     WRITE(6,*)I,FPTR(0),FPTR(1),FPTR(2)
  ENDDO
ENDIF

! find area of each box side that intersects with the tetrahedron

CALL GET_AREAS(VERTS,FACES,NFACES,WHICH_SIDE,AREAS)
RETURN
END SUBROUTINE GET_VERTS

!  ------------------ GET_VERTS2D ------------------------

SUBROUTINE GET_VERTS2D(XB,PLANE_INDEX,VV0,VV1,VV2,VERTS,NVERTS)

! find vertices formed by the intersection of a rectangle and a triangle.
! The rectangle and triangle lie in the same plane.  The rectangle sides
! are aligned with the coordinate axes
!                        2
!                       . .
!                     .   .
!   2---------------.--3  .
!   |              .   |   .
!   |             0    |   .
!   |                . |   .
!   |                  | .  .
!   |                  |    1
!   0------------------1

REAL(EB), DIMENSION(0:5), INTENT(IN) :: XB ! xmin, xmax, ymin, ymax, zmin, zmax
INTEGER, INTENT(IN) :: PLANE_INDEX         ! 0, 1, or 2
REAL(EB), DIMENSION(0:2), INTENT(IN) :: VV0, VV1, VV2
REAL(EB), INTENT(OUT), TARGET :: VERTS(0:299)
INTEGER, INTENT(OUT) :: NVERTS

REAL(EB), DIMENSION(0:4) :: RECT
REAL(EB), TARGET :: TRI_VERTS(0:5), RECT_VERTS(0:7)
REAL(EB), POINTER, DIMENSION(:) :: V0, V1, V2
INTEGER, TARGET :: TRI_EDGES(0:5)=(/0,1, 1,2,  2,0/)
INTEGER, TARGET :: RECT_EDGES(0:7)=(/0,1, 1,3, 3,2, 2,0 /)
INTEGER, POINTER, DIMENSION(:) :: RE, TE

REAL(EB) :: PLANE_VAL
INTEGER :: I, J
REAL(EB), DIMENSION(0:199), TARGET :: VERTS2D
REAL(EB), POINTER, DIMENSION(:) :: V, V2D
REAL(EB), POINTER, DIMENSION(:) :: VIN, VOUT
REAL(EB), POINTER, DIMENSION(:) :: R0, R1, T0, T1
REAL(EB), TARGET :: VERT_SEGS(0:3)
INTEGER :: NVERT_SEGS
REAL(EB) :: VMID(0:1), VCOPY(0:1)

NVERTS=0
PLANE_VAL = XB(PLANE_INDEX)
IF (PLANE_INDEX==0) THEN
  TRI_VERTS(0:1) = VV0(1:2)
  TRI_VERTS(2:3) = VV1(1:2)
  TRI_VERTS(4:5) = VV2(1:2)
  RECT(0:3) = XB(2:5)
ELSE IF (PLANE_INDEX==1) THEN
  TRI_VERTS(0:1) = VV0(1:2)
  TRI_VERTS(2:3) = VV1(1:2)
  TRI_VERTS(4:5) = VV2(1:2)
  RECT(0:3) = XB(2:5)
ELSE IF (PLANE_INDEX==2) THEN
  TRI_VERTS(0:1) = (/VV0(0),VV0(2)/)
  TRI_VERTS(2:3) = (/VV1(0),VV1(2)/)
  TRI_VERTS(4:5) = (/VV2(0),VV2(2)/)
  RECT(0:3) = (/XB(0:1),XB(4:5)/)
ELSE IF (PLANE_INDEX==3) THEN
  TRI_VERTS(0:1) = (/VV0(0),VV0(2)/)
  TRI_VERTS(2:3) = (/VV1(0),VV1(2)/)
  TRI_VERTS(4:5) = (/VV2(0),VV2(2)/)
  RECT(0:3) = (/XB(0:1),XB(4:5)/)
ELSE IF (PLANE_INDEX==4) THEN
  TRI_VERTS(0:1) = VV0(0:1)
  TRI_VERTS(2:3) = VV1(0:1)
  TRI_VERTS(4:5) = VV2(0:1)
  RECT(0:3) = XB(0:3)
ELSE
  TRI_VERTS(0:1) = VV0(0:1)
  TRI_VERTS(2:3) = VV1(0:1)
  TRI_VERTS(4:5) = VV2(0:1)
  RECT(0:3) = XB(0:3)
ENDIF
RECT_VERTS(0:1)=(/RECT(0),RECT(2)/)
RECT_VERTS(2:3)=(/RECT(1),RECT(2)/)
RECT_VERTS(4:5)=(/RECT(0),RECT(3)/)
RECT_VERTS(6:7)=(/RECT(1),RECT(3)/)
V0(0:1)=>TRI_VERTS(0:1)
V1(0:1)=>TRI_VERTS(2:3)
V2(0:1)=>TRI_VERTS(4:5)
VMID=(V0+V1+V2)/3.0_EB
IF (.NOT.IN_TRIANGLE2D(V0,V1,V2,VMID)) THEN ! make sure vertices are ordered counter-clockwise
   VCOPY=TRI_VERTS(0:1)
   TRI_VERTS(0:1)=TRI_VERTS(2:3)
   TRI_VERTS(2:3)=VCOPY
   V0(0:1)=>TRI_VERTS(0:1)
   V1(0:1)=>TRI_VERTS(2:3)
ENDIF

! check for triangle verts inside rectangles

DO I = 0, 2
   VIN(0:1)=>TRI_VERTS(2*I:2*I+1)
   IF (IN_RECTANGLE2D(RECT,VIN)) THEN
      VOUT(0:1)=>VERTS2D(2*NVERTS:2*NVERTS+1)
      VOUT(0:1) = VIN(0:1)
      NVERTS = NVERTS + 1
   ENDIF
ENDDO

! check for rectangle verts inside triangles

DO I = 0, 3
   VIN(0:1)=>RECT_VERTS(2*I:2*I+1)
   IF (IN_TRIANGLE2D(V0,V1,V2,VIN)) THEN
      VOUT(0:1)=>VERTS2D(2*NVERTS:2*NVERTS+1)
      VOUT(0:1) = VIN(0:1)
      NVERTS = NVERTS + 1
   ENDIF
ENDDO

! check for rectangle edges that intersect with triangle edges

DO I = 0, 3 ! rectangle edges
   RE(0:1)=>RECT_EDGES(2*I:2*I+1)
   R0(0:1)=>RECT_VERTS(2*RE(0):2*RE(0)+1)
   R1(0:1)=>RECT_VERTS(2*RE(1):2*RE(1)+1)
   DO J = 0, 2 ! triangle edges
      TE(0:1)=>TRI_EDGES(2*J:2*J+1)
      T0(0:1)=>TRI_VERTS(2*TE(0):2*TE(0)+1)
      T1(0:1)=>TRI_VERTS(2*TE(1):2*TE(1)+1)
      IF (LINE_SEGMENT_INTERSECT(R0,R1,T0,T1,VERT_SEGS,NVERT_SEGS)) THEN
         VIN(0:1)=>VERT_SEGS(0:1)
         VOUT(0:1)=>VERTS2D(2*NVERTS:2*NVERTS+1)
         VOUT(0:1)=VIN(0:1)
         NVERTS=NVERTS+1

         IF (NVERT_SEGS>1) THEN
            VIN(0:1)=>VERT_SEGS(2:3)
            VOUT(0:1)=>VERTS2D(2*NVERTS:2*NVERTS+1)
            VOUT(0:1)=VIN(0:1)
            NVERTS=NVERTS+1
         ENDIF
      ENDIF
   ENDDO
ENDDO

! copy 2d vertex info to 3d vertex array

DO I = 0, NVERTS-1
   V(0:2)=>VERTS(3*I:3*I+2)
   V2D(0:1)=>VERTS2D(2*I:2*I+1)
   IF (PLANE_INDEX==0 .OR. PLANE_INDEX==1) THEN
      V(0:2)=(/PLANE_VAL,V2D(0:1)/)
   ELSE IF (PLANE_INDEX==2 .OR. PLANE_INDEX==3) THEN
      V(0:2)=(/V2D(0),PLANE_VAL,V2D(1)/)
   ELSE
      V(0:2)=(/V2D(0:1),PLANE_VAL/)
   ENDIF
ENDDO

END SUBROUTINE GET_VERTS2D

!  ------------------ REMOVE_DEGENERATE_FACES ------------------------

SUBROUTINE REMOVE_DEGENERATE_FACES(N_FACES,MAX_FACES,FACES,WHICH_SIDE)
INTEGER, INTENT(INOUT) :: N_FACES
INTEGER, INTENT(IN) :: MAX_FACES
INTEGER, TARGET, INTENT(INOUT) :: FACES(3*MAX_FACES), WHICH_SIDE(MAX_FACES)

INTEGER, POINTER, DIMENSION(:) :: FACEFROM, FACETO
INTEGER :: IFROM, N_FACES2

IFROM = 1
N_FACES2=0
DO IFROM = 1, N_FACES
   FACEFROM(1:3) => FACES(3*IFROM-2:3*IFROM)
   IF (FACEFROM(1)==FACEFROM(2) .OR. FACEFROM(1)==FACEFROM(3) .OR. FACEFROM(2)==FACEFROM(3)) CYCLE
   N_FACES2=N_FACES2+1
   IF (N_FACES2==IFROM) CYCLE
   FACETO(1:3) => FACES(3*N_FACES2-2:3*N_FACES2)
   FACETO(1:3) = FACEFROM(1:3)
   WHICH_SIDE(N_FACES2) = WHICH_SIDE(IFROM)
ENDDO
N_FACES=N_FACES2

END SUBROUTINE REMOVE_DEGENERATE_FACES

! ---------------------------- REMOVE_DUPLICATE_VERTS ----------------------------------------

SUBROUTINE REMOVE_DUPLICATE_VERTS(N_VERTS,N_FACES,N_VOLUS,&
                            MAX_VERTS,MAX_FACES,MAX_VOLUS,FIRST_FACE_INDEX,&
                            VERTS,FACES,VOLUS,VERT_EPS)
INTEGER, INTENT(INOUT) :: N_VERTS
INTEGER, INTENT(IN) :: N_FACES, N_VOLUS
INTEGER, INTENT(IN) :: MAX_VERTS, MAX_FACES, MAX_VOLUS
INTEGER, INTENT(IN) :: FIRST_FACE_INDEX
INTEGER, TARGET, INTENT(INOUT) :: VOLUS(4*MAX_VOLUS)
REAL(EB), TARGET, INTENT(INOUT) :: VERTS(3*MAX_VERTS)
INTEGER, TARGET, INTENT(INOUT) :: FACES(3*MAX_FACES)
REAL(EB), INTENT(IN) :: VERT_EPS

REAL(EB), POINTER, DIMENSION(:) :: VI, VJ
REAL(EB), DIMENSION(3) :: VIMVJ
REAL(EB) :: NORM_VI, NORM_VJ, NORM_VIMVJ

INTEGER :: I, J, K
INTEGER :: OFFSET

OFFSET = 1-FIRST_FACE_INDEX
I = 1
DO WHILE (I<=N_VERTS)
   VI=>VERTS(3*I-2:3*I)
   NORM_VI = MAX(ABS(VI(1)),ABS(VI(2)),ABS(VI(3)))
   J = I+1
   DO WHILE (J<=N_VERTS)
      VJ=>VERTS(3*J-2:3*J)
      VIMVJ = VI-VJ
      NORM_VJ = MAX(ABS(VJ(1)),ABS(VJ(2)),ABS(VJ(3)))
      NORM_VIMVJ = MAX(ABS(VIMVJ(1)),ABS(VIMVJ(2)),ABS(VIMVJ(3)))
      IF (NORM_VIMVJ <= MAX(1.0_EB,NORM_VI,NORM_VJ)*VERT_EPS) THEN
         ! vertex I and J are the same
         ! first copy index J -> I in face list
         ! next copy index N_VERTS -> J in face list
         ! finally reduce N_VERTS by 1

         DO K = 1, 3*N_FACES
           IF (FACES(K)==J-OFFSET)FACES(K)=I-OFFSET
         ENDDO
         DO K = 1, 4*N_VOLUS
           IF (VOLUS(K)==J-OFFSET)VOLUS(K)=I-OFFSET
         ENDDO
         VJ(1:3)=VERTS(3*N_VERTS-2:3*N_VERTS)

         DO K = 1, 3*N_FACES
           IF (FACES(K)==N_VERTS-OFFSET)FACES(K)=J-OFFSET
         ENDDO
         DO K = 1, 4*N_VOLUS
           IF (VOLUS(K)==N_VERTS-OFFSET)VOLUS(K)=J-OFFSET
         ENDDO
         N_VERTS=N_VERTS-1
         CYCLE
      ENDIF
      J=J+1
   ENDDO
   I=I+1
ENDDO
END SUBROUTINE REMOVE_DUPLICATE_VERTS

!  ------------------ IN_RECTANGLE2D ------------------------

LOGICAL FUNCTION LINE_SEGMENT_INTERSECT(R0,R1,T0,T1,VERT_SEGS,NVERT_SEGS)
REAL(EB), INTENT(IN), DIMENSION(0:1) :: R0, R1, T0, T1
REAL(EB), INTENT(OUT) :: VERT_SEGS(0:3)
INTEGER, INTENT(OUT) :: NVERT_SEGS
REAL(EB) :: A(2,2), B(2), U(2), DENOM
REAL(EB), PARAMETER :: EPS=0.0001_EB
REAL(EB) :: RXMIN, RXMAX, RYMIN, RYMAX
REAL(EB) :: TXMIN, TXMAX, TYMIN, TYMAX

NVERT_SEGS=0
LINE_SEGMENT_INTERSECT = .FALSE.

! handle case when segments coincide vertically
! (there are two solutions in this case)

IF (ABS(T0(0)-T1(0))<EPS .AND. ABS(R0(0)-R1(0))<EPS) THEN
   TYMIN = MIN(T0(1),T1(1))
   TYMAX = MIN(T0(1),T1(1))
   RYMIN = MIN(R0(1),R1(1))
   RYMAX = MIN(R0(1),R1(1))
   IF (MIN(TYMAX,RYMAX)>MAX(TYMIN,RYMIN)) THEN
      VERT_SEGS(0)=T0(0)
      VERT_SEGS(1)=MIN(TYMAX,RYMAX)

      VERT_SEGS(2)=T0(0)
      VERT_SEGS(3)=MAX(TYMIN,RYMIN)
      LINE_SEGMENT_INTERSECT = .TRUE.
      NVERT_SEGS=2
   ENDIF
   RETURN
ENDIF

! handle case when segments coincide horizontally
! (there are two solutions in this case)

IF (ABS(T0(1)-T1(1))<EPS .AND. ABS(R0(1)-R1(1))<EPS) THEN
   TXMIN = MIN(T0(0),T1(0))
   TXMAX = MIN(T0(0),T1(0))
   RXMIN = MIN(R0(0),R1(0))
   RXMAX = MIN(R0(0),R1(0))
   IF (MIN(TXMAX,RXMAX)>MAX(TXMIN,RXMIN)) THEN
      VERT_SEGS(0)=MIN(TXMAX,RXMAX)
      VERT_SEGS(1)=T0(1)

      VERT_SEGS(2)=MAX(TXMIN,RXMIN)
      VERT_SEGS(3)=T0(1)
      LINE_SEGMENT_INTERSECT = .TRUE.
      NVERT_SEGS=2
   ENDIF
   RETURN
ENDIF

! find where two lines cross by solving a 2x2 system of linear equations
! the intersection occurs within the two line segments if and only if
! the solution U(1), U(2) lies between 0 and 1

A(1,1) = R1(0)-R0(0)
A(1,2) = -(T1(0)-T0(0))
A(2,1) = R1(1)-R0(1)!
A(2,2) = -(T1(1)-T0(1))

B(1) = (T0(0)-R0(0))
B(2) = (T0(1)-R0(1))

DENOM = A(1,1)*A(2,2)-A(1,2)*A(2,1)
IF (ABS(DENOM)<EPS) RETURN
U(1)=(B(1)*A(2,2)-B(2)*A(1,2))/DENOM
U(2)=(A(1,1)*B(2)-A(2,1)*B(1))/DENOM

! only consider solution when 0<=U<=1

IF (U(1)<0.0_EB .OR. U(1)>1.0_EB .OR. U(2)<0.0_EB .OR. U(2)>1.0_EB) RETURN

VERT_SEGS(0)=R0(0)+(R1(0)-R0(0))*U(1)
VERT_SEGS(1)=R0(1)+(R1(1)-R0(1))*U(1)
NVERT_SEGS=1
LINE_SEGMENT_INTERSECT = .TRUE.
END FUNCTION LINE_SEGMENT_INTERSECT

!  ------------------ IN_RECTANGLE2D ------------------------

LOGICAL FUNCTION IN_RECTANGLE2D(XB,VERT)
! tests whether VERT is inside 2d rectangle defined by XB
REAL(EB), INTENT(IN) :: XB(0:3)
REAL(EB), INTENT(IN) :: VERT(0:1)

IF (VERT(0)<XB(0) .OR. VERT(0)>XB(1) .OR. &
   VERT(1)<XB(2) .OR. VERT(1)>XB(3)) THEN
   IN_RECTANGLE2D=.FALSE.
ELSE
   IN_RECTANGLE2D=.TRUE.
ENDIF

END FUNCTION IN_RECTANGLE2D

!  ------------------ IN_TRIANGLE2D ------------------------

LOGICAL FUNCTION IN_TRIANGLE2D(V0,V1,V2,VERT)
REAL(EB), DIMENSION(0:1), INTENT(IN) :: V0, V1, V2, VERT

REAL(EB), DIMENSION(0:1) :: DV, DVERT, N

!        V0
!       / \
!      /   \
!     /     \
!    /       \
!   /         \
!  V1---------V2

IN_TRIANGLE2D=.FALSE.

DV=V1-V0
DVERT=VERT-V0
N(0:1)=(/DV(1),-DV(0)/)
IF (N(0)*DVERT(0)+N(1)*DVERT(1)>0.0_EB) RETURN

DV=V2-V1
DVERT=VERT-V1
N(0:1)=(/DV(1),-DV(0)/)
IF (N(0)*DVERT(0)+N(1)*DVERT(1)>0.0_EB) RETURN

DV=V0-V2
DVERT=VERT-V2
N(0:1)=(/DV(1),-DV(0)/)
IF (N(0)*DVERT(0)+N(1)*DVERT(1)>0.0_EB) RETURN

IN_TRIANGLE2D=.TRUE.
END FUNCTION IN_TRIANGLE2D

!  ------------------ COMPARE_VERTS ------------------------

INTEGER FUNCTION COMPARE_VERTS(A,B,C,N)

! determine orientation of vertices A, B, C relative to vector N

REAL(EB), DIMENSION(0:2), INTENT(IN) :: A, B, C, N
REAL(EB), DIMENSION(0:2) :: AMC, BMC, ACROSSB
REAL(EB) :: VAL

AMC = A - C
BMC = B - C
CALL CROSS_PRODUCT(ACROSSB,AMC, BMC)
VAL = DOT_PRODUCT(ACROSSB,N)/6.0_EB
IF (VAL>0.0_EB) THEN
  COMPARE_VERTS = 1
ELSE IF (VAL==0.0_EB) THEN
  COMPARE_VERTS = 0
ELSE
  COMPARE_VERTS = -1
ENDIF
END FUNCTION COMPARE_VERTS

!  ------------------ ORDER_VERTS ------------------------

SUBROUTINE ORDER_VERTS(FACEVERTS,NVERTS,SIDE_INDEX)

! order vertices of a given face

INTEGER, INTENT(IN) :: NVERTS, SIDE_INDEX
REAL(EB), DIMENSION(0:3*NVERTS-1), INTENT(INOUT), TARGET :: FACEVERTS

REAL(EB), DIMENSION(:), POINTER :: NORMAL
REAL(EB), DIMENSION(0:2) :: VERT_CENTER
REAL(EB), DIMENSION(:), POINTER :: VERT, VERT1, VERT2
REAL(EB), DIMENSION(0:3*NVERTS-1) :: FACEVERTS_TEMP
INTEGER :: ORIENTATION
INTEGER :: I, II, IIP1, J
INTEGER, DIMENSION(0:100) :: ORDER
INTEGER :: N_CLOCKWISE, N_COUNTER_CLOCKWISE

! first 6 (starting at 0) normal indices are for box faces
! next 4 (6->9) normal indices are for tetrahedron faces

IF (SIDE_INDEX<=5) THEN
   NORMAL(0:2) => BOX_NORMALS(0:2,SIDE_INDEX)
ELSE
   NORMAL(0:2) => TETRA_NORMALS(0:2,SIDE_INDEX-6)
ENDIF

! find center of vertices

VERT_CENTER=0.0_EB
DO I = 0, NVERTS-1
  VERT(0:2) => FACEVERTS(3*I:3*I+2)
  VERT_CENTER = VERT_CENTER + VERT
ENDDO
VERT_CENTER = VERT_CENTER/REAL(NVERTS,EB)

! split vertices into two parts
!    those before and those after vertex 0

N_CLOCKWISE=0
N_COUNTER_CLOCKWISE=0
VERT1(0:2) => FACEVERTS(0:2)
DO I = 1, NVERTS-1
  VERT2(0:2) => FACEVERTS(3*I:3*I+2)
  ORIENTATION=COMPARE_VERTS(VERT1,VERT2,VERT_CENTER,NORMAL)
  IF (ORIENTATION>=1) THEN
    ORDER(NVERTS-1-N_CLOCKWISE) = I
    N_CLOCKWISE = N_CLOCKWISE + 1
  ELSE
    ORDER(N_COUNTER_CLOCKWISE) = I
    N_COUNTER_CLOCKWISE = N_COUNTER_CLOCKWISE + 1
  ENDIF
ENDDO
ORDER(N_COUNTER_CLOCKWISE) = 0

! order vertices that are counter clockwise from vertex 0

IF (N_COUNTER_CLOCKWISE>1) THEN
   DO J = 0, N_COUNTER_CLOCKWISE - 1
   DO I = 0, N_COUNTER_CLOCKWISE - 2
      II = ORDER(I)
      VERT1(0:2) => FACEVERTS(3*II:3*II+2)
      IIP1 = ORDER(I+1)
      VERT2(0:2) => FACEVERTS(3*IIP1:3*IIP1+2)
      ORIENTATION=COMPARE_VERTS(VERT1,VERT2,VERT_CENTER,NORMAL)
      IF (ORIENTATION<1) THEN
         ORDER(I) = IIP1
         ORDER(I+1) = II
      ENDIF
   ENDDO
   ENDDO
ENDIF

! order vertices that are counter clockwise from vertex 0

IF (N_CLOCKWISE>1) THEN
   DO J = 0, N_CLOCKWISE - 1
   DO I = 0, N_CLOCKWISE - 2
      II = ORDER(N_COUNTER_CLOCKWISE+1+I)
      VERT1(0:2) => FACEVERTS(3*II:3*II+2)
      IIP1 = ORDER(N_COUNTER_CLOCKWISE+1+I+1)
      VERT2(0:2) => FACEVERTS(3*IIP1:3*IIP1+2)
      ORIENTATION=COMPARE_VERTS(VERT1,VERT2,VERT_CENTER,NORMAL)
      IF (ORIENTATION<1) THEN
         ORDER(N_COUNTER_CLOCKWISE+1+I) = IIP1
         ORDER(N_COUNTER_CLOCKWISE+1+I+1) = II
      ENDIF
   ENDDO
   ENDDO
ENDIF

! copy ordered vertices into original data structure

FACEVERTS_TEMP(0:3*NVERTS-1) = FACEVERTS(0:3*NVERTS-1)
DO I = 0, NVERTS-1
   FACEVERTS(3*I:3*I+2) = FACEVERTS_TEMP(3*ORDER(I):3*ORDER(I)+2)
ENDDO

END SUBROUTINE ORDER_VERTS

!  ------------------ SETUP_VERTS ------------------------

SUBROUTINE SETUP_VERTS(BOX_BOUNDS,V0,V1,V2,V3)
REAL(EB), DIMENSION(0:5), INTENT(IN) :: BOX_BOUNDS
REAL(EB), DIMENSION(0:2), INTENT(IN) :: V0, V1, V2, V3
REAL(EB), DIMENSION(0:2) :: VDIFF1, VDIFF2, VCROSS
INTEGER :: TP, BP
INTEGER, DIMENSION(:), POINTER :: VERTS

! define box vertices given x/y/z min/max values

BOX_VERTS(0,0) = BOX_BOUNDS(MIN_X)
BOX_VERTS(1,0) = BOX_BOUNDS(MIN_Y)
BOX_VERTS(2,0) = BOX_BOUNDS(MIN_Z)

BOX_VERTS(0,1) = BOX_BOUNDS(MAX_X)
BOX_VERTS(1,1) = BOX_BOUNDS(MIN_Y)
BOX_VERTS(2,1) = BOX_BOUNDS(MIN_Z)

BOX_VERTS(0,2) = BOX_BOUNDS(MIN_X)
BOX_VERTS(1,2) = BOX_BOUNDS(MAX_Y)
BOX_VERTS(2,2) = BOX_BOUNDS(MIN_Z)

BOX_VERTS(0,3) = BOX_BOUNDS(MAX_X)
BOX_VERTS(1,3) = BOX_BOUNDS(MAX_Y)
BOX_VERTS(2,3) = BOX_BOUNDS(MIN_Z)

BOX_VERTS(0,4) = BOX_BOUNDS(MIN_X)
BOX_VERTS(1,4) = BOX_BOUNDS(MIN_Y)
BOX_VERTS(2,4) = BOX_BOUNDS(MAX_Z)

BOX_VERTS(0,5) = BOX_BOUNDS(MAX_X)
BOX_VERTS(1,5) = BOX_BOUNDS(MIN_Y)
BOX_VERTS(2,5) = BOX_BOUNDS(MAX_Z)

BOX_VERTS(0,6) = BOX_BOUNDS(MIN_X)
BOX_VERTS(1,6) = BOX_BOUNDS(MAX_Y)
BOX_VERTS(2,6) = BOX_BOUNDS(MAX_Z)

BOX_VERTS(0,7) = BOX_BOUNDS(MAX_X)
BOX_VERTS(1,7) = BOX_BOUNDS(MAX_Y)
BOX_VERTS(2,7) = BOX_BOUNDS(MAX_Z)

TETRA_BOUNDS(0) = MIN(V0(0),V1(0),V2(0),V3(0))
TETRA_BOUNDS(1) = MAX(V0(0),V1(0),V2(0),V3(0))
TETRA_BOUNDS(2) = MIN(V0(1),V1(1),V2(1),V3(1))
TETRA_BOUNDS(3) = MAX(V0(1),V1(1),V2(1),V3(1))
TETRA_BOUNDS(4) = MIN(V0(2),V1(2),V2(2),V3(2))
TETRA_BOUNDS(5) = MAX(V0(2),V1(2),V2(2),V3(2))

DO TP = 0, 3
  VERTS(0:2) => TETRA_PLANE2VERT(0:2,TP)
  TETRA_PLANE_BOUNDS(0,TP) = MIN(TETRA_VERTS(0,VERTS(0)),TETRA_VERTS(0,VERTS(1)),TETRA_VERTS(0,VERTS(2)))
  TETRA_PLANE_BOUNDS(1,TP) = MAX(TETRA_VERTS(0,VERTS(0)),TETRA_VERTS(0,VERTS(1)),TETRA_VERTS(0,VERTS(2)))
  TETRA_PLANE_BOUNDS(2,TP) = MIN(TETRA_VERTS(1,VERTS(0)),TETRA_VERTS(1,VERTS(1)),TETRA_VERTS(1,VERTS(2)))
  TETRA_PLANE_BOUNDS(3,TP) = MAX(TETRA_VERTS(1,VERTS(0)),TETRA_VERTS(1,VERTS(1)),TETRA_VERTS(1,VERTS(2)))
  TETRA_PLANE_BOUNDS(4,TP) = MIN(TETRA_VERTS(2,VERTS(0)),TETRA_VERTS(2,VERTS(1)),TETRA_VERTS(2,VERTS(2)))
  TETRA_PLANE_BOUNDS(5,TP) = MAX(TETRA_VERTS(2,VERTS(0)),TETRA_VERTS(2,VERTS(1)),TETRA_VERTS(2,VERTS(2)))
ENDDO

DO BP = 0, 5
  VERTS(0:3) => BOX_PLANE2VERT(0:3,BP)
  BOX_PLANE_BOUNDS(0,BP) = MIN(BOX_VERTS(0,VERTS(0)),BOX_VERTS(0,VERTS(1)),BOX_VERTS(0,VERTS(2)),BOX_VERTS(0,VERTS(3)))
  BOX_PLANE_BOUNDS(1,BP) = MAX(BOX_VERTS(0,VERTS(0)),BOX_VERTS(0,VERTS(1)),BOX_VERTS(0,VERTS(2)),BOX_VERTS(0,VERTS(3)))
  BOX_PLANE_BOUNDS(2,BP) = MIN(BOX_VERTS(1,VERTS(0)),BOX_VERTS(1,VERTS(1)),BOX_VERTS(1,VERTS(2)),BOX_VERTS(0,VERTS(3)))
  BOX_PLANE_BOUNDS(3,BP) = MAX(BOX_VERTS(1,VERTS(0)),BOX_VERTS(1,VERTS(1)),BOX_VERTS(1,VERTS(2)),BOX_VERTS(0,VERTS(3)))
  BOX_PLANE_BOUNDS(4,BP) = MIN(BOX_VERTS(2,VERTS(0)),BOX_VERTS(2,VERTS(1)),BOX_VERTS(2,VERTS(2)),BOX_VERTS(0,VERTS(3)))
  BOX_PLANE_BOUNDS(5,BP) = MAX(BOX_VERTS(2,VERTS(0)),BOX_VERTS(2,VERTS(1)),BOX_VERTS(2,VERTS(2)),BOX_VERTS(0,VERTS(3)))
ENDDO

! define tetrahedron vertices

TETRA_VERTS(0:2,0) = V0
TETRA_VERTS(0:2,1) = V1
TETRA_VERTS(0:2,2) = V2
TETRA_VERTS(0:2,3) = V3

!           3
!          /.\
!         / . \
!        /  .  \
!       /   .   \
!      /    2    \
!     /   .   .   \
!    /  .       .  \
!   / .           . \
!  0-----------------1

! compute normals for each tetrahedron face
! orient vertices (clockwise) so that normals point outwards

VDIFF1 = V1 - V0
VDIFF2 = V3 - V0
CALL CROSS_PRODUCT(VCROSS,VDIFF1,VDIFF2)
TETRA_NORMALS(0:2,0)=VEC_NORMALIZE(VCROSS)

VDIFF1 = V2 - V1
VDIFF2 = V3 - V1
CALL CROSS_PRODUCT(VCROSS,VDIFF1,VDIFF2)
TETRA_NORMALS(0:2,1)=VEC_NORMALIZE(VCROSS)

VDIFF1 = V0 - V2
VDIFF2 = V3 - V2
CALL CROSS_PRODUCT(VCROSS,VDIFF1,VDIFF2)
TETRA_NORMALS(0:2,2)=VEC_NORMALIZE(VCROSS)

VDIFF1 = V0 - V1
VDIFF2 = V2 - V1
CALL CROSS_PRODUCT(VCROSS,VDIFF1,VDIFF2)
TETRA_NORMALS(0:2,3)=VEC_NORMALIZE(VCROSS)

RETURN
END SUBROUTINE SETUP_VERTS

!  ------------------ PLANE_EDGE_INTERSECTION ------------------------

INTEGER FUNCTION PLANE_EDGE_INTERSECTION(X0,N0,V0,V1,VERT)

! PLANE: (X-X0).DOT.N0 = 0, EDGE: V0->V1

! find T such that ( (1-T)*V0 + T*V1 - X0 ) .DOT. N0 = 0
! if V1-V0 is perpendicular to N0 ie (V1-V0).DOT.N0 = 0 then there is no solution

REAL(EB), DIMENSION(0:2), INTENT(IN) :: X0, N0, V0, V1
REAL(EB), DIMENSION(0:2), INTENT(OUT) :: VERT
REAL(EB), DIMENSION(0:2) :: V1MV0, V0MX0

REAL(EB) :: DENOM, T

PLANE_EDGE_INTERSECTION=0

V1MV0 = V1-V0
V0MX0 = V0-X0
DENOM = DOT_PRODUCT(V1MV0,N0)
IF (DENOM/=0) THEN
   T = -DOT_PRODUCT(V0MX0,N0)/DENOM
   IF (T>=0.0_EB .AND.T<=1.0_EB) THEN
      PLANE_EDGE_INTERSECTION=1
      VERT = V0 + T*V1MV0
   ENDIF
ENDIF
RETURN
END FUNCTION PLANE_EDGE_INTERSECTION

!  ------------------ BOXPLANE_EDGE_INTERSECTION ------------------------

INTEGER FUNCTION BOXPLANE_EDGE_INTERSECTION(X0,PLANE,V0,V1,VERT)

! PLANE: (X-X0).DOT.N0 = 0, EDGE: V0->V1

! find T such that ( (1-T)*V0 + T*V1 - X0 ) .DOT. N0 = 0
! if V1-V0 is perpendicular to N0 ie (V1-V0).DOT.N0 = 0 then there is no solution

REAL(EB), DIMENSION(0:2), INTENT(IN) :: X0, V0, V1
INTEGER, INTENT(IN) :: PLANE
REAL(EB), DIMENSION(0:2), INTENT(OUT) :: VERT
REAL(EB), DIMENSION(0:2) :: V1MV0, V0MX0

REAL(EB) :: DENOM, T
INTEGER, DIMENSION(0:5) :: PLANE_INDEX

DATA PLANE_INDEX/0,0,1,1,2,2/

BOXPLANE_EDGE_INTERSECTION=0

V1MV0 = V1-V0
V0MX0 = V0-X0
!DENOM = DOT_PRODUCT(V1MV0,N0)
DENOM = V1MV0(PLANE_INDEX(PLANE))
IF (DENOM/=0) THEN
!   T = -DOT_PRODUCT(V0MX0,N0)/DENOM
   T = -V0MX0(PLANE_INDEX(PLANE))/DENOM
   IF (T>=0.0_EB .AND.T<=1.0_EB) THEN
      BOXPLANE_EDGE_INTERSECTION=1
      VERT = V0 + T*V1MV0
   ENDIF
ENDIF
RETURN
END FUNCTION BOXPLANE_EDGE_INTERSECTION

!  ------------------ IN_BOX ------------------------

INTEGER FUNCTION IN_BOX(XYZ, IGNORE_PLANE)

! determine if a vertex XYZ is inside the box

REAL(EB), INTENT(IN), DIMENSION(0:2) :: XYZ
INTEGER, INTENT(IN) :: IGNORE_PLANE

INTEGER BP
REAL(EB), DIMENSION(0:2) :: VECDIFF
REAL(EB), DIMENSION(:), POINTER :: BOXVERT, BOXNORMAL

IN_BOX=1
DO BP = 0, 5
   IF (BP==IGNORE_PLANE) CYCLE
   BOXVERT(0:2) => BOX_VERTS(0:2,BOX_PLANE2VERT(0,BP))
   BOXNORMAL(0:2) => BOX_NORMALS(0:2,BP)
   VECDIFF = XYZ - BOXVERT
   IF (DOT_PRODUCT(BOXNORMAL,VECDIFF)>0.0_EB) THEN
      IN_BOX=0
      RETURN
   ENDIF
ENDDO
END FUNCTION IN_BOX

!  ------------------ IN_TETRA ------------------------

INTEGER FUNCTION IN_TETRA(XYZ, IGNORE_PLANE)

! determine if a vertex XYZ is inside the tetrahedron

REAL(EB), INTENT(IN), DIMENSION(0:2) :: XYZ
INTEGER, INTENT(IN) :: IGNORE_PLANE
REAL(EB), DIMENSION(0:2) :: VECDIFF
REAL(EB), DIMENSION(:), POINTER :: TETRAVERT, TETRANORMAL

INTEGER :: TP

IN_TETRA=1
DO TP = 0, 3
   IF (TP==IGNORE_PLANE) CYCLE
   TETRAVERT(0:2) => TETRA_VERTS(0:2,TETRA_PLANE2VERT(0,TP))
   TETRANORMAL(0:2) => TETRA_NORMALS(0:2,TP)
   VECDIFF = XYZ - TETRAVERT
   IF (DOT_PRODUCT(TETRANORMAL,VECDIFF)>0.0_EB) THEN
      IN_TETRA=0
      RETURN
   ENDIF
ENDDO
END FUNCTION IN_TETRA

!  ------------------ TEST_IN_TETRA0 ------------------------

SUBROUTINE TEST_IN_TETRA0(XYZ, IN_TETRA, TETRA_STATE)

! determine if a vertex XYZ is inside the tetrahedron

REAL(EB), INTENT(IN), DIMENSION(0:2) :: XYZ
INTEGER, INTENT(OUT) :: TETRA_STATE(0:3), IN_TETRA

REAL(EB), DIMENSION(0:2) :: VECDIFF
REAL(EB), DIMENSION(:), POINTER :: TETRAVERT, TETRANORMAL

INTEGER :: TP

IN_TETRA=1
DO TP = 0, 3
   TETRAVERT(0:2) => TETRA_VERTS(0:2,TETRA_PLANE2VERT(0,TP))
   TETRANORMAL(0:2) => TETRA_NORMALS(0:2,TP)
   VECDIFF = XYZ - TETRAVERT
   IF (DOT_PRODUCT(TETRANORMAL,VECDIFF)>0.0_EB) THEN
      TETRA_STATE(TP) = 0
      IN_TETRA=0
   ELSE
      TETRA_STATE(TP) = 1
   ENDIF
ENDDO
END SUBROUTINE TEST_IN_TETRA0

!  ------------------ VEC_NORMALIZE ------------------------

FUNCTION VEC_NORMALIZE(U)

! normalize a vector so |U|=1

REAL(EB), DIMENSION(0:2) :: VEC_NORMALIZE
REAL(EB), DIMENSION(0:2), INTENT(IN) :: U

REAL(EB) :: SUM

SUM = SQRT(DOT_PRODUCT(U,U))
IF (SUM/=0.0_EB) THEN
  VEC_NORMALIZE = U/SUM
ELSE
  VEC_NORMALIZE = U
ENDIF

RETURN
END FUNCTION VEC_NORMALIZE

!  ------------------ ON_BOX_PLANE ------------------------

INTEGER FUNCTION ON_BOX_PLANE(XB,V1,V2,V3)
REAL(EB), INTENT(IN) :: XB(6), V1(3), V2(3), V3(3)
INTEGER, PARAMETER :: XMIN=1,XMAX=2,YMIN=3,YMAX=4,ZMIN=6,ZMAX=6
INTEGER, PARAMETER :: X=1, Y=2, Z=3
REAL(EB), PARAMETER :: EPS_ABS=0.001_EB
REAL(EB) :: EPS

! determines whether the plane formed by the vertices V1, V2 and V3 coincides with one
! of the 6 box sides defined by XB

EPS = EPS_ABS*MAX(1.0_EB,XB(XMIN))
IF (ABS(V1(X)-XB(XMIN))<EPS .AND. ABS(V2(X)-XB(XMIN))<EPS .AND. ABS(V3(X)-XB(XMIN))<EPS) THEN
   ON_BOX_PLANE=0
   RETURN
ENDIF

EPS = EPS_ABS*MAX(1.0_EB,XB(XMAX))
IF (ABS(V1(X)-XB(XMAX))<EPS .AND. ABS(V2(X)-XB(XMAX))<EPS .AND. ABS(V3(X)-XB(XMAX))<EPS) THEN
   ON_BOX_PLANE=1
   RETURN
ENDIF

EPS = EPS_ABS*MAX(1.0_EB,XB(YMIN))
IF (ABS(V1(Y)-XB(YMIN))<EPS .AND. ABS(V2(Y)-XB(YMIN))<EPS .AND. ABS(V3(Y)-XB(YMIN))<EPS) THEN
   ON_BOX_PLANE=2
   RETURN
ENDIF

EPS = EPS_ABS*MAX(1.0_EB,XB(YMAX))
IF (ABS(V1(Y)-XB(YMAX))<EPS .AND. ABS(V2(Y)-XB(YMAX))<EPS .AND. ABS(V3(Y)-XB(YMAX))<EPS) THEN
   ON_BOX_PLANE=3
   RETURN
ENDIF

EPS = EPS_ABS*MAX(1.0_EB,XB(ZMIN))
IF (ABS(V1(Z)-XB(ZMIN))<EPS .AND. ABS(V2(Z)-XB(ZMIN))<EPS .AND. ABS(V3(Z)-XB(ZMIN))<EPS) THEN
   ON_BOX_PLANE=4
   RETURN
ENDIF

EPS = EPS_ABS*MAX(1.0_EB,XB(ZMAX))
IF (ABS(V1(Z)-XB(ZMAX))<EPS .AND. ABS(V2(Z)-XB(ZMAX))<EPS .AND. ABS(V3(Z)-XB(ZMAX))<EPS) THEN
   ON_BOX_PLANE=5
   RETURN
ENDIF
ON_BOX_PLANE=-1
END FUNCTION ON_BOX_PLANE

#ifdef MISC_GSMV_GEOM
!  ------------------ RANK3 ------------------------

SUBROUTINE RANK3(V1,V2,V3,I1,I2,I3)
INTEGER, INTENT(IN) :: V1, V2, V3
INTEGER, INTENT(OUT) :: I1, I2, I3

I1 = MIN(V1,V2,V3)
I3 = MAX(V1,V2,V3)
I2 = V1+V2+V3-I1-I3

END SUBROUTINE RANK3

!  ------------------ VOLS2FACES ------------------------

SUBROUTINE VOLS2SURFACE(VOLS,NVOLS,FACES,NFACES)

! find the exterior surface formed from a collection of volumes

INTEGER, INTENT(IN) :: NVOLS
INTEGER, INTENT(IN) :: VOLS(NVOLS,4)
INTEGER, INTENT(OUT) :: NFACES
INTEGER, INTENT(OUT), DIMENSION(4*NVOLS,3) :: FACES

INTEGER :: I, ML, MM, MU
INTEGER :: FACEDUP(4*NVOLS)

INTEGER :: ITO, IFROM

NFACES = 4*NVOLS
DO I = 1, NVOLS
   CALL RANK3(VOLS(I,4),VOLS(I,1),VOLS(I,2),ML,MM,MU)
   FACES(4*I-3,1) = ML
   FACES(4*I-3,2) = MM
   FACES(4*I-3,3) = MU

   CALL RANK3(VOLS(I,4),VOLS(I,2),VOLS(I,3),ML,MM,MU)
   FACES(4*I-2,1) = ML
   FACES(4*I-2,2) = MM
   FACES(4*I-2,3) = MU

   CALL RANK3(VOLS(I,4),VOLS(I,3),VOLS(I,1),ML,MM,MU)
   FACES(4*I-1,1) = ML
   FACES(4*I-1,2) = MM
   FACES(4*I-1,3) = MU

   CALL RANK3(VOLS(I,1),VOLS(I,3),VOLS(I,2),ML,MM,MU)
   FACES(4*I  ,1) = ML
   FACES(4*I  ,2) = MM
   FACES(4*I  ,3) = MU
ENDDO

! tag interior faces (faces connected to more than one volume)

FACEDUP(1:NFACES)=0
DO I = 1, NFACES
   IF (FACEDUP(I)==1) CYCLE
   DO J = I+1, NFACES
      IF (FACEDUP(J)==1) CYCLE
      IF (FACES(I,1)==FACES(J,1) .AND.FACES(I,2)==FACES(J,2) .AND.FACES(I,3)==FACES(J,3))FACEDUP(J)=1
   ENDDO
ENDDO

! remove duplicate (ie interior) faces

ITO=1
DO IFROM = 2, NFACES
   IF (FACEDUP(IFROM)==1) CYCLE
   ITO = ITO + 1
   IF (ITO/=I) THEN
     FACES(ITO,1) = FACES(IFROM,1)
     FACES(ITO,2) = FACES(IFROM,2)
     FACES(ITO,3) = FACES(IFROM,3)
   ENDIF
ENDDO
NFACES=ITO

END SUBROUTINE VOLS2SURFACE

!  ------------------ FACES2EDGES ------------------------

SUBROUTINE FACES2EDGES(FACES,NFACES,EDGES,NEDGES)

! determine edges from faces, eliminate all duplicate edges

INTEGER, INTENT(IN) :: NFACES
INTEGER, INTENT(IN), DIMENSION(NFACES,3) :: FACES
INTEGER, INTENT(OUT), DIMENSION(3*NFACES,2) :: EDGES
INTEGER, INTENT(OUT) :: NEDGES

INTEGER :: I, TEMP1, TEMP2, FINISHED, IFROM, ITO

NEDGES=3*NFACES
DO I = 1, NFACES
   EDGES(3*I-2,1) = MIN(FACES(I,1),FACES(I,2))
   EDGES(3*I-2,2) = MAX(FACES(I,1),FACES(I,2))
   EDGES(3*I-1,1) = MIN(FACES(I,1),FACES(I,3))
   EDGES(3*I-1,2) = MAX(FACES(I,1),FACES(I,3))
   EDGES(3*I  ,1) = MIN(FACES(I,2),FACES(I,3))
   EDGES(3*I  ,2) = MAX(FACES(I,2),FACES(I,3))
ENDDO

! this sort is o(n**2), replace with an o(n*log(n)) sort once n is "large"

FINISHED=1
DO WHILE (FINISHED==1)
   DO I = 2, NEDGES
      IF (EDGES(I-1,1)<EDGES(I,1)) CYCLE
      IF (EDGES(I-1,1)==EDGES(I,1) .AND.EDGES(I-1,2)<EDGES(I,2)) CYCLE
      TEMP1=EDGES(I-1,1)
      TEMP2=EDGES(I-1,2)
      EDGES(I-1,1)=EDGES(I,1)
      EDGES(I-1,2)=EDGES(I,2)
      EDGES(I,1)=TEMP1
      EDGES(I,2)=TEMP2
      FINISHED=0
   ENDDO
ENDDO

! remove duplicates

ITO=1
DO IFROM = 2, NEDGES
   IF (EDGES(IFROM,1)==EDGES(IFROM-1,1) .AND.EDGES(IFROM,1)==EDGES(IFROM-1,2)) CYCLE
   ITO = ITO + 1
   IF (ITO/=I) THEN
     EDGES(ITO,1) = EDGES(IFROM,1)
     EDGES(ITO,2) = EDGES(IFROM,2)
   ENDIF
ENDDO
NEDGES = ITO

END SUBROUTINE FACES2EDGES

!  ------------------ CLASSIFY_GEOM ------------------------

SUBROUTINE CLASSIFY_GEOM(XGRID,NX,YGRID,NY,ZGRID,NZ,FACES,NFACES,VERTS,NVERTS,GEOM_STATE)

! identify which grid cells are cut cells
! a grid cell is a cut cell if at least on triangle in FACES intersects with the grid cell

INTEGER, INTENT(IN) :: NX, NY, NZ
REAL(EB), INTENT(IN) :: XGRID(NX), YGRID(NY), ZGRID(NZ)
INTEGER, INTENT(IN) :: NFACES, NVERTS
INTEGER, INTENT(IN), DIMENSION(NFACES) :: FACES(NFACES,3)
REAL(EB), INTENT(IN), DIMENSION(NVERTS,3) :: VERTS
INTEGER, INTENT(OUT), DIMENSION(NX-1,NY-1,NZ-1) :: GEOM_STATE

INTEGER :: I,J,K,L
REAL(EB) :: XB(6),V1(3),V2(3),V3(3)
INTEGER :: STATE

DO I = 1, NX-1
   XB(1:2) = XGRID(I:I+1)
   DO J = 1, NY-1
     XB(3:4) = YGRID(J:J+1)
      DO K = 1, NZ-1
         XB(5:6) = ZGRID(K:K+1)
         DO L = 1, NFACES
            V1(1:3)=VERTS(FACES(I,1),1:3)
            V2(1:3)=VERTS(FACES(I,2),1:3)
            V3(1:3)=VERTS(FACES(I,3),1:3)
            STATE = BOXTRIANGLE(XB,V1,V2,V3)
            IF (STATE==1)EXIT
         ENDDO
         GEOM_STATE(I,J,K)=STATE
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE CLASSIFY_GEOM
#endif

!  ------------------ IN_BOX3 ------------------------

INTEGER FUNCTION IN_BOX3(XB,V)
REAL(EB) :: XB(6), V(3)

IF (XB(1)<=V(1) .AND.V(1)<=XB(2) .AND.&
    XB(3)<=V(2) .AND.V(2)<=XB(4) .AND.&
    XB(5)<=V(3) .AND.V(3)<=XB(6)) THEN
   IN_BOX3=1
   RETURN
ENDIF
IN_BOX3=0
END FUNCTION

!  ------------------ IN_BOX2 ------------------------

INTEGER FUNCTION IN_BOX2(XB,V)
REAL(EB) :: XB(6), V(3)

IF (XB(1)<V(1) .AND.V(1)<XB(2) .AND.&
   XB(3)<V(2) .AND.V(2)<XB(4) .AND.&
   XB(5)<V(3) .AND.V(3)<XB(6)) THEN
   IN_BOX2=1
   RETURN
ENDIF
IN_BOX2=0
END FUNCTION

!  ------------------ BOXTRIANGLE ------------------------

INTEGER FUNCTION BOXTRIANGLE(XB,V1,V2,V3)

! determine if a box XB intersects with a triangle (V1, V2, V3)

REAL(EB), INTENT(IN) :: XB(6), V1(3), V2(3), V3(3)

IF (ON_BOX_PLANE(XB,V1,V2,V3)==-1) THEN
   BOXTRIANGLE=0
   RETURN
ENDIF
IF (IN_BOX2(XB,V1)==1 .OR. IN_BOX2(XB,V2)==1 .OR. IN_BOX2(XB,V3)==1) THEN
   BOXTRIANGLE=1
   RETURN
ENDIF

BOXTRIANGLE=0

END FUNCTION BOXTRIANGLE

!  ------------------ DISTANCE3 ------------------------

REAL(EB) FUNCTION DISTANCE3(V1,V2)
REAL(EB), INTENT(IN), DIMENSION(3) :: V1, V2
REAL(EB) :: DX, DY, DZ

DX = V1(1)-V2(1)
DY = V1(2)-V2(2)
DZ = V1(3)-V2(3)
DISTANCE3 = SQRT(DX*DX+DY*DY+DZ*DZ)
END FUNCTION DISTANCE3

!  ------------------ GET_VERTTYPE ------------------------

SUBROUTINE GET_VERTTYPE(NVERTS, TRIANGLES, NTRIANGLES, VERT_TYPE)
INTEGER, INTENT(IN) :: NVERTS, NTRIANGLES
INTEGER, INTENT(OUT), DIMENSION(NVERTS) :: VERT_TYPE
INTEGER, INTENT(IN), DIMENSION(3*NTRIANGLES), TARGET :: TRIANGLES

! classify each vertex in a geometry as either interior or exterior
! a vertex VI is interior if the vertices connected to I form a cycle or loop
! ie VI is connected to v1, v2, v3 and v1 -> v2 -> v3 -> v1
! if they don't form a loop then it is an exterior vertex

! exterior vertices (connected to a blockage or mesh boundary) won't be moved or deleted

INTEGER, DIMENSION(NVERTS) :: VERT_COUNT, TRIANGLE_COUNT
INTEGER :: I, J, K, VERTJ_INDEX, TRIJ_INDEX, MAXCOUNT, VERTK_INDEX
INTEGER, POINTER, DIMENSION(:) :: TRII
INTEGER, DIMENSION(:,:), ALLOCATABLE :: VERT_TRILIST

! count number of triangles connected to each vertex

TRIANGLE_COUNT(1:NVERTS) = 0
DO I = 1, NTRIANGLES
   TRII=>TRIANGLES(3*I-2:3*I)
   DO J = 1, 3
     VERTJ_INDEX = TRII(J)
     IF(VERTJ_INDEX>=1 .AND. VERTJ_INDEX<=NVERTS)TRIANGLE_COUNT(VERTJ_INDEX) = TRIANGLE_COUNT(VERTJ_INDEX) + 1
   END DO
END DO

MAXCOUNT = TRIANGLE_COUNT(1)
DO I = 2, NVERTS
  MAXCOUNT = MAX(MAXCOUNT,TRIANGLE_COUNT(I))
END DO

! construct a list of triangles connected to each vertex
! VERT_TRILIST(I,1) contains number of triangles connected to vertex I
! VERT_TRILIST(I,2-> ) contains the triangle indices

ALLOCATE(VERT_TRILIST(NVERTS,MAXCOUNT+1))
VERT_TRILIST(1:NVERTS,1)=0
DO I = 1, NTRIANGLES
   TRII=>TRIANGLES(3*I-2:3*I)
   DO J = 1, 3
     VERTJ_INDEX = TRII(J)
     IF(VERTJ_INDEX>=1 .AND. VERTJ_INDEX<=NVERTS)THEN
        VERT_TRILIST(VERTJ_INDEX,1) = VERT_TRILIST(VERTJ_INDEX,1) + 1 ! bump triangle count by 1 for vertex j
        K = VERT_TRILIST(VERTJ_INDEX,1)+1
        VERT_TRILIST(VERTJ_INDEX,K) = I  ! put triangle index into triangle count + 1
     ENDIF
   END DO
END DO

! count vertices connected to each vertex

VERT_TYPE(1:NVERTS) = 1
OUTER: DO I = 1, NVERTS
  VERT_COUNT(1:NVERTS) = 0
  DO J = 1, VERT_TRILIST(I,1) ! loop over triangles connected to vertex I
    TRIJ_INDEX = VERT_TRILIST(I,1+J)
    DO K = 1, 3 ! loop over vertices of triangle J
      VERTK_INDEX = TRIANGLES(3*TRIJ_INDEX-3+K)
      IF(VERTK_INDEX/=I)VERT_COUNT(VERTK_INDEX) = VERT_COUNT(VERTK_INDEX) + 1
    END DO
  END DO
  DO J = 1, NVERTS
    IF(VERT_COUNT(J)==1)THEN ! consider all vertices that are connected to vertex I 
      VERT_TYPE(I) = 0       ! if all of these neighbors have two neighbors among this set then I is interior
      CYCLE OUTER            ! if at least one of these neighbors has only one neigbor among this set then I is on the exterior
    ENDIF
  END DO
END DO OUTER

DEALLOCATE(VERT_TRILIST)

END SUBROUTINE GET_VERTTYPE

!  ------------------ DECIMATE_FB ------------------------

SUBROUTINE DECIMATE_FB(VERTS, NVERTS, FACES, NFACES, MESH_BOUNDS_FB, DELTA)
INTEGER, INTENT(INOUT) :: NVERTS, NFACES
REAL(FB), INTENT(INOUT), DIMENSION(3*NVERTS) :: VERTS
INTEGER, INTENT(INOUT), DIMENSION(3*NFACES) :: FACES
REAL(FB), INTENT(IN) :: DELTA
REAL(FB), DIMENSION(6), INTENT(IN) :: MESH_BOUNDS_FB

REAL(EB), DIMENSION(3*NVERTS) :: VERTS_EB
REAL(EB) :: DELTA_EB
REAL(EB), DIMENSION(6) :: MESH_BOUNDS_EB

DELTA_EB = REAL(DELTA,EB)
VERTS_EB(1:3*NVERTS) = REAL(VERTS(1:3*NVERTS),EB)
MESH_BOUNDS_EB(1:6) = REAL(MESH_BOUNDS_FB(1:6),EB)
CALL DECIMATE(VERTS_EB,NVERTS,FACES,NFACES,MESH_BOUNDS_EB,DELTA_EB)
VERTS(1:3*NVERTS) = REAL(VERTS_EB(1:3*NVERTS),FB)
RETURN
END SUBROUTINE DECIMATE_FB

!  ------------------ AVERAGE_VERTS2 ------------------------

SUBROUTINE AVERAGE_VERTS2(V1,V1TYPE,V2,V2TYPE,MESH_BOUNDS,VAVG)
REAL(EB), DIMENSION(3), INTENT(IN) :: V1, V2
INTEGER, INTENT(IN) :: V1TYPE, V2TYPE
REAL(EB), DIMENSION(6), INTENT(IN) :: MESH_BOUNDS
REAL(EB), DIMENSION(3), INTENT(OUT) :: VAVG

REAL(EB) :: BOXEPS

BOXEPS = 0.001

IF (V1TYPE == 0) THEN
  VAVG(1:3) = V1(1:3)
  RETURN
ENDIF

IF (V2TYPE == 0) THEN
  VAVG(1:3) = V2(1:3)
  RETURN
ENDIF

IF (ABS(V1(1)-MESH_BOUNDS(1))<BOXEPS .OR. ABS(V1(1)-MESH_BOUNDS(2))<BOXEPS) THEN
  VAVG(1) = V1(1)
ELSE IF (ABS(V2(1)-MESH_BOUNDS(1))<BOXEPS .OR. ABS(V2(1)-MESH_BOUNDS(2))<BOXEPS) THEN
  VAVG(1) = V2(1)
ELSE
  VAVG(1) = (V1(1)+V2(1))/2.0
ENDIF

IF (ABS(V1(2)-MESH_BOUNDS(3))<BOXEPS .OR. ABS(V1(2)-MESH_BOUNDS(4))<BOXEPS) THEN
  VAVG(2) = V1(2)
ELSE IF (ABS(V2(2)-MESH_BOUNDS(3))<BOXEPS .OR. ABS(V2(2)-MESH_BOUNDS(4))<BOXEPS) THEN
  VAVG(2) = V2(2)
ELSE
  VAVG(2) = (V1(2)+V2(2))/2.0
ENDIF

IF ( ABS( V1(3)-MESH_BOUNDS(5) )<BOXEPS .OR. ABS(V1(3)-MESH_BOUNDS(6))<BOXEPS) THEN
  VAVG(3) = V1(3)
ELSE IF (ABS(V2(3)-MESH_BOUNDS(5))<BOXEPS .OR. ABS(V2(3)-MESH_BOUNDS(6))<BOXEPS) THEN
  VAVG(3) = V2(3)
ELSE
  VAVG(3) = (V1(3)+V2(3))/2.0
ENDIF

END SUBROUTINE AVERAGE_VERTS2

!  ------------------ AVERAGE_VERTS3 ------------------------

SUBROUTINE AVERAGE_VERTS3(V1,V1TYPE,V2,V2TYPE,V3,V3TYPE,MESH_BOUNDS,VAVG)
REAL(EB), DIMENSION(3), INTENT(IN) :: V1, V2, V3
INTEGER, INTENT(IN) :: V1TYPE, V2TYPE, V3TYPE
REAL(EB), DIMENSION(6), INTENT(IN) :: MESH_BOUNDS
REAL(EB), DIMENSION(3), INTENT(OUT) :: VAVG

REAL(EB) :: BOXEPS

BOXEPS = 0.001

IF (V1TYPE==0) THEN
  VAVG(1:3) = V1(1:3)
  RETURN
ENDIF

IF (V2TYPE==0) THEN
  VAVG(1:3) = V2(1:3)
  RETURN
ENDIF

IF (V3TYPE==0) THEN
  VAVG(1:3) = V3(1:3)
  RETURN
ENDIF

IF (ABS(V1(1)-MESH_BOUNDS(1))<BOXEPS .OR. ABS(V1(1)-MESH_BOUNDS(2))<BOXEPS) THEN
  VAVG(1) = V1(1)
ELSE IF (ABS(V2(1)-MESH_BOUNDS(1))<BOXEPS .OR. ABS(V2(1)-MESH_BOUNDS(2))<BOXEPS) THEN
  VAVG(1) = V2(1)
ELSE IF (ABS(V3(1)-MESH_BOUNDS(1))<BOXEPS .OR. ABS(V3(1)-MESH_BOUNDS(2))<BOXEPS) THEN
  VAVG(1) = V3(1)
ELSE
  VAVG(1) = (V1(1)+V2(1)+V3(1))/3.0
ENDIF

IF (ABS(V1(2)-MESH_BOUNDS(3))<BOXEPS .OR. ABS(V1(2)-MESH_BOUNDS(4))<BOXEPS) THEN
  VAVG(2) = V1(2)
ELSE IF (ABS(V2(2)-MESH_BOUNDS(3))<BOXEPS .OR. ABS(V2(2)-MESH_BOUNDS(4))<BOXEPS) THEN
  VAVG(2) = V2(2)
ELSE IF (ABS(V3(2)-MESH_BOUNDS(3))<BOXEPS .OR. ABS(V3(2)-MESH_BOUNDS(4))<BOXEPS) THEN
  VAVG(2) = V3(2)
ELSE
  VAVG(2) = (V1(2)+V2(2)+V3(2))/3.0
ENDIF

IF (ABS(V1(3)-MESH_BOUNDS(5))<BOXEPS .OR. ABS(V1(3)-MESH_BOUNDS(6))<BOXEPS) THEN
  VAVG(3) = V1(3)
ELSE IF (ABS(V2(3)-MESH_BOUNDS(5))<BOXEPS .OR. ABS(V2(3)-MESH_BOUNDS(6))<BOXEPS) THEN
  VAVG(3) = V2(3)
ELSE IF (ABS(V3(3)-MESH_BOUNDS(5))<BOXEPS .OR. ABS(V3(3)-MESH_BOUNDS(6))<BOXEPS) THEN
  VAVG(3) = V3(3)
ELSE
  VAVG(3) = (V1(3)+V2(3)+V3(3))/3.0
ENDIF

END SUBROUTINE AVERAGE_VERTS3

!  ------------------ DECIMATE ------------------------

SUBROUTINE DECIMATE(VERTS, NVERTS, FACES, NFACES, MESH_BOUNDS, DELTA)
! This routine reduces the size of a geometry by
!  1) merging vertices that are "close" together
!  2) eliminating redundent vertices
!  3) eliminating "singular" triangles

INTEGER, INTENT(INOUT) :: NVERTS, NFACES
REAL(EB), INTENT(INOUT), DIMENSION(3*NVERTS), TARGET :: VERTS
INTEGER, INTENT(INOUT), DIMENSION(3*NFACES), TARGET :: FACES
REAL(EB), INTENT(IN) :: DELTA
REAL(EB), DIMENSION(6), INTENT(IN) :: MESH_BOUNDS

INTEGER, PARAMETER :: V_MERGED=-1, V_DISCARD=0, V_ORIGINAL=1

INTEGER, DIMENSION(NVERTS) :: VERT_STATE, VERT_MAP, VERT_TYPE
REAL(EB), POINTER, DIMENSION(:) :: V1, V2, V3, VERTFROM, VERTTO
INTEGER, POINTER, DIMENSION(:) :: TRI_I, TRI_FROM, TRI_TO
INTEGER, DIMENSION(3) :: TRI_NEW
REAL(EB) :: D12, D13, D23
INTEGER :: I, IFROM, ITO, ITER, MAX_ITER
LOGICAL :: HAVE_SMALL
REAL(EB), DIMENSION(3) :: VAVG

HAVE_SMALL = .TRUE.
MAX_ITER = 4
ITER = 0
DO WHILE (HAVE_SMALL .AND. ITER<MAX_ITER) ! iterate until no further changes are made (or 10 times whichever comes first)
   HAVE_SMALL = .FALSE.

   ! VERT_STATE
   !    V_MERGE =   -1  -  merged vertex
   !    V_DISCARD =  0  -  discard vertex
   !    V_ORIGINAL = 1  -  vertex kept and not changed

   VERT_STATE(1:NVERTS) = V_ORIGINAL

   DO I = 1, NVERTS
     VERT_MAP(I) = I
   END DO

   ITER = ITER + 1

   CALL GET_VERTTYPE(NVERTS, FACES, NFACES, VERT_TYPE)

! combine vertices that are close together

   TRI_LOOP: DO I = 1, NFACES
      TRI_I=>FACES(3*I-2:3*I)
      V1=>VERTS(3*TRI_I(1)-2:3*TRI_I(1))
      V2=>VERTS(3*TRI_I(2)-2:3*TRI_I(2))
      V3=>VERTS(3*TRI_I(3)-2:3*TRI_I(3))

      IF (VERT_STATE(TRI_I(1))/=V_ORIGINAL .OR. &  ! only look at triangles that have not changed
          VERT_STATE(TRI_I(2))/=V_ORIGINAL .OR. &
          VERT_STATE(TRI_I(3))/=V_ORIGINAL) CYCLE TRI_LOOP

      D12 = DISTANCE3(V1,V2)
      D13 = DISTANCE3(V1,V3)
      D23 = DISTANCE3(V2,V3)
      IF (D12>DELTA .AND.D13>DELTA .AND.D23>DELTA) CYCLE TRI_LOOP ! triangle too large, do not combine verts

      HAVE_SMALL = .TRUE.
      IF (D12<DELTA .AND.D13>DELTA .AND.D23>DELTA) THEN ! combine verts 1 and 2 leave 3 alone
         VERT_STATE(TRI_I(1)) = V_MERGED
         VERT_STATE(TRI_I(2)) = V_DISCARD
         CALL AVERAGE_VERTS2(V1,VERT_TYPE(TRI_I(1)),V2,VERT_TYPE(TRI_I(2)),MESH_BOUNDS,VAVG)
         V1(1:3) = VAVG(1:3)
         VERT_MAP(TRI_I(2)) = TRI_I(1)
         TRI_I(2) = TRI_I(1)
      ELSE IF (D13<DELTA .AND.D12>DELTA .AND.D23>DELTA) THEN ! combine verts 1 and 3
         VERT_STATE(TRI_I(1)) = V_MERGED
         VERT_STATE(TRI_I(3)) = V_DISCARD
         CALL AVERAGE_VERTS2(V1,VERT_TYPE(TRI_I(1)),V3,VERT_TYPE(TRI_I(3)),MESH_BOUNDS,VAVG)
         V1(1:3) = VAVG(1:3)
         VERT_MAP(TRI_I(3)) = TRI_I(1)
         TRI_I(3) = TRI_I(1)
      ELSE IF (D23<DELTA .AND.D12>DELTA .AND.D13>DELTA) THEN ! combine verts 2 and 3
         VERT_STATE(TRI_I(2)) = V_MERGED
         VERT_STATE(TRI_I(3)) = V_DISCARD
         CALL AVERAGE_VERTS2(V2,VERT_TYPE(TRI_I(2)),V3,VERT_TYPE(TRI_I(3)),MESH_BOUNDS,VAVG)
         V2(1:3) = VAVG(1:3)
         VERT_MAP(TRI_I(3)) = TRI_I(2)
         TRI_I(3) = TRI_I(2)
      ELSE  ! combine verts 1, 2 and 3
         VERT_STATE(TRI_I(1))= V_MERGED
         VERT_STATE(TRI_I(2))= V_DISCARD
         VERT_STATE(TRI_I(3))= V_DISCARD
         CALL AVERAGE_VERTS3(V1, VERT_TYPE(TRI_I(1)), V2, VERT_TYPE(TRI_I(2)), V3, VERT_TYPE(TRI_I(3)), MESH_BOUNDS,VAVG)
         V1(1:3) = VAVG(1:3)
         VERT_MAP(TRI_I(2)) = TRI_I(1)
         VERT_MAP(TRI_I(3)) = TRI_I(1)
         TRI_I(2) = TRI_I(1)
         TRI_I(3) = TRI_I(1)
      ENDIF
   ENDDO TRI_LOOP

   ! remap triangle vertices

   DO I = 1, NFACES
      TRI_I=>FACES(3*I-2:3*I)
      TRI_I(1) = VERT_MAP(TRI_I(1))
      TRI_I(2) = VERT_MAP(TRI_I(2))
      TRI_I(3) = VERT_MAP(TRI_I(3))
   END DO

   ! construct new vertex list skipping over vertices that have been removed

   ITO = 0
   DO IFROM = 1, NVERTS
      IF (VERT_STATE(IFROM)/=V_DISCARD) THEN
        ITO = ITO + 1
        VERTFROM=>VERTS(3*IFROM-2:3*IFROM)
        VERTTO=>VERTS(3*ITO-2:3*ITO)
        VERTTO(1:3)=VERTFROM(1:3)
      ENDIF
      VERT_MAP(IFROM) = ITO
   ENDDO
   NVERTS = ITO

   ! eliminate singular triangles (as a result of merged vertices)

   ITO = 0
   DO IFROM = 1, NFACES
      TRI_FROM=>FACES(3*IFROM-2:3*IFROM)
      IF (TRI_FROM(1)/=TRI_FROM(2) .AND. TRI_FROM(1)/=TRI_FROM(3) .AND. TRI_FROM(2)/=TRI_FROM(3)) THEN
        TRI_NEW(1)=VERT_MAP(TRI_FROM(1))
        TRI_NEW(2)=VERT_MAP(TRI_FROM(2))
        TRI_NEW(3)=VERT_MAP(TRI_FROM(3))
        IF (TRI_NEW(1)/=0 .AND. TRI_NEW(2)/=0 .AND. TRI_NEW(3)/=0) THEN
          ITO=ITO+1
          TRI_TO=>FACES(3*ITO-2:3*ITO)
          TRI_TO(1)=TRI_NEW(1)
          TRI_TO(2)=TRI_NEW(2)
          TRI_TO(3)=TRI_NEW(3)
        ENDIF
      ENDIF
   ENDDO
   NFACES = ITO
ENDDO

END SUBROUTINE DECIMATE

!  ------------------ GEOMCLIP ------------------------

SUBROUTINE GEOMCLIP(VERTS, NVERTS, FACES, NFACES, XB)
INTEGER, INTENT(INOUT) :: NVERTS, NFACES
REAL(EB), INTENT(INOUT), DIMENSION(3*NVERTS), TARGET :: VERTS
INTEGER, INTENT(INOUT), DIMENSION(3*NFACES), TARGET :: FACES
REAL(EB), DIMENSION(6), INTENT(IN) :: XB

INTEGER, POINTER, DIMENSION(:) :: TRII
REAL(EB), POINTER, DIMENSION(:) :: V1, V2, V3
INTEGER, DIMENSION(NVERTS) :: HAVE_VERT, VERT_MAP
INTEGER, DIMENSION(NFACES) :: HAVE_TRIANGLE
INTEGER :: NT, NV

! find triangles that are not within the box defined by XB

HAVE_VERT(1:NVERTS) = 0
HAVE_TRIANGLE(1:NFACES) = 0
DO I = 1, NFACES
  TRII(1:3) => FACES(3*I-2:3*I)
  V1(1:3) => VERTS(3*TRII(1)-2:3*TRII(1))
  V2(1:3) => VERTS(3*TRII(2)-2:3*TRII(2))
  V3(1:3) => VERTS(3*TRII(3)-2:3*TRII(3))
  IF (IN_BOX3(XB,V1)==1 .OR. IN_BOX3(XB,V2)==1 .OR. IN_BOX3(XB,V3)==1) THEN
    HAVE_VERT(TRII(1)) = 1
    HAVE_VERT(TRII(2)) = 1
    HAVE_VERT(TRII(3)) = 1
    HAVE_TRIANGLE(I) = 1
  ENDIF
END DO

! construct a map pointing from old vertex index to new vertex index

NV = 0
DO I = 1, NVERTS
  IF (HAVE_VERT(I) == 1) THEN
    NV = NV + 1
    IF (I /= NV) THEN
       VERTS(3*NV-2:3*NV) = VERTS(3*I-2:3*I)
    ENDIF
  ENDIF
  VERT_MAP(I) = NV
END DO

! define new vertex indices

NT=0
DO I = 1, NFACES
  IF (HAVE_TRIANGLE(I)==1) THEN
    TRII(1:3) => FACES(3*I-2:3*I)
    TRII(1) = VERT_MAP(TRII(1))
    TRII(2) = VERT_MAP(TRII(2))
    TRII(3) = VERT_MAP(TRII(3))
    NT = NT + 1
    IF (NT /= I)FACES(3*NT-2:3*NT) = TRII(1:3)
  ENDIF
END DO

NFACES = NT
NVERTS = NV

END SUBROUTINE GEOMCLIP

!  ------------------ IN_TRIANGLE ------------------------

LOGICAL FUNCTION IN_TRIANGLE(VERT,V1,V2,V3,ON_PERIMETER)

! determine whether VERT is inside the triangle defined by V1, V2, V3 (in a plane)

REAL(EB), INTENT(IN), DIMENSION(2) :: VERT, V1, V2, V3
INTEGER, INTENT(IN) :: ON_PERIMETER

REAL(EB), DIMENSION(2) :: VH, V2H, V3H
REAL(EB) :: DET, A, B

! reject if VERT is outside the bounding box for V1, V2 and V3

IN_TRIANGLE = .FALSE.
IF ( VERT(1)<MIN(V1(1),V2(1),V3(1)) .OR. VERT(2)<MIN(V1(2),V2(2),V3(2)).OR.&
     VERT(1)>MAX(V1(1),V2(1),V3(1)) .OR. VERT(2)>MAX(V1(2),V2(2),V3(2))) RETURN

! VERT = V1 + a(V2-V1) + b(V3-V1)

! rewrite as
! ( V2H_x   V3H_x )  (a)     (  VH_x )
! (               )  ( )  =  (       )
! ( V2H_y   VH3_y )  (b)     (  VH_y )

! where
! VH = VERT - V1
! V2H = V2 - V1
! V3H = V3 - V1

! solution:
!     | VH_x   V3H_x |     |  V2H_x  V3H_x |
! a = |              |  /  |               |
!     | VH_y   VH3_y |     |  V2H_y  VH3_y |

!     | V2H_x  VH_x  |     |  V2H_x  V3H_x |
! b = |              |  /  |               |
!     | V2H_y  VH_y  |     |  V2H_y  VH3_y |


! VERT is in triangle if 0<=a<=1 and 0<=b<=1 and 0<=a+b<=1

 VH = VERT - V1
V2H = V2 - V1
V3H = V3 - V1


DET = V2H(1)*V3H(2) - V2H(2)*V3H(1)
IF ( DET/=0.0_EB) THEN
   A = ( VH(1)*V3H(2) -  VH(2)*V3H(1))/DET
   B = (V2H(1)* VH(2) - V2H(2)* VH(1))/DET

   IF (ON_PERIMETER == 0) THEN
! reject if VERT is outside of triangle or on its perimeter
      IF (A<=0.0_EB.OR.B<=0.0_EB.OR.A+B>=1.0_EB) RETURN
   ELSE
! reject if VERT is outside of triangle
      IF (A<0.0_EB.OR.B<0.0_EB.OR.A+B>1.0_EB) RETURN
   ENDIF
ENDIF
IN_TRIANGLE = .TRUE.
END FUNCTION

!  ------------------ IN_POLYGON ------------------------

LOGICAL FUNCTION IN_POLYGON(VERT,VERTS,NVERTS,TRIANGLES,NTRIANGLES)
REAL(EB), INTENT(IN) :: VERT(1:3)
INTEGER, INTENT(IN) :: NVERTS
REAL(EB), INTENT(IN), DIMENSION(3*NVERTS), TARGET :: VERTS
INTEGER, INTENT(IN) :: NTRIANGLES
INTEGER, INTENT(IN) :: TRIANGLES(3*NTRIANGLES)

INTEGER :: I, I1, I2, I3
REAL(EB), DIMENSION(:), POINTER :: V1, V2, V3

IN_POLYGON = .TRUE.
DO I = 1, NTRIANGLES
    I1 = TRIANGLES(3*I-2)
    I2 = TRIANGLES(3*I-1)
    I3 = TRIANGLES(3*I)
    V1 => VERTS(3*I1-2:3*I1)
    V2 => VERTS(3*I2-2:3*I2)
    V3 => VERTS(3*I3-2:3*I3)
    IF(IN_TRIANGLE(VERT,V1,V2,V3,1))RETURN
END DO
IN_POLYGON = .FALSE.
END FUNCTION IN_POLYGON

!  ------------------ IS_ACUTE_ANGLE ------------------------

LOGICAL FUNCTION IS_ANGLE_GE_180(V1,V2,V3)

! determine whether the angle formed by vertices V1, V2, V3 is acute ( <180.0) or obtuse (>= 180.0)

REAL(EB), INTENT(IN), DIMENSION(2) :: V1, V2, V3
REAL(EB), DIMENSION(2) :: DV1, DV2, DV1PERP

IS_ANGLE_GE_180 = .TRUE.
DV1 = V2 - V1
DV1PERP = (/-DV1(2),DV1(1)/)
DV2 = V3 - V2
IF ( DOT_PRODUCT(DV1PERP,DV2) > 0.0_EB) IS_ANGLE_GE_180 = .FALSE.
END FUNCTION IS_ANGLE_GE_180

!  ------------------ POLY2TRI ------------------------

SUBROUTINE POLY2TRI(VERTS,NVERTS,POLY,NPOLY,TRIS,NTRIS)

! given a list of vertices forming a polygon, create a 2nd list of triangles

INTEGER, INTENT(IN) :: NVERTS, NPOLY
REAL(EB), INTENT(IN), DIMENSION(2*NVERTS), TARGET :: VERTS
INTEGER, INTENT(IN), DIMENSION(NPOLY) :: POLY

INTEGER, INTENT(OUT) :: NTRIS
INTEGER, INTENT(OUT), DIMENSION(3*(NPOLY-2)) :: TRIS

REAL(EB), POINTER, DIMENSION(:) :: VERT, V1, V2, V3
INTEGER, DIMENSION(NPOLY) :: P
INTEGER :: NP, I, J, IP

NP = NPOLY
P(1:NP) = POLY(1:NPOLY)
NTRIS = 0
IP = 1
DO WHILE (NP >= 3 .AND. IP .LE. NPOLY)
   LOOP_I: DO I = 1, NP - 2

      ! candidate triangle
      V1 => VERTS(2*P(I  )-1:2*P(I  ))
      V2 => VERTS(2*P(I+1)-1:2*P(I+1))
      V3 => VERTS(2*P(I+2)-1:2*P(I+2))

      IF(IS_ANGLE_GE_180(V1,V2,V3))CYCLE LOOP_I ! reject triangle if angle >= 180.0
      DO J = 1, I-1
         VERT => VERTS(2*POLY(J)-1:2*POLY(J))        ! reject triangle if another polygon point is inside this triangle
         IF (IN_TRIANGLE(VERT,V1,V2,V3,0)) CYCLE LOOP_I
      END DO
      ! skip over verts in candidate triangle, ie I, I+1 and I+2
      DO J = I+3, NP
         VERT => VERTS(2*POLY(J)-1:2*POLY(J))        ! reject triangle if another polygon point is inside this triangle
         IF (IN_TRIANGLE(VERT,V1,V2,V3,0)) CYCLE LOOP_I
      END DO
      ! V1, V2, V3 is a valid triangle so copy it into the triangle list
      TRIS(3*NTRIS+1) = P(I)
      TRIS(3*NTRIS+2) = P(I+1)
      TRIS(3*NTRIS+3) = P(I+2)
      NTRIS = NTRIS + 1

      ! reduce polygon count by 1 (take out polygon pt corresponding to V2 )
      DO J = I + 2, NP
         P(J-1) = P(J)
      END DO
      NP = NP - 1
      EXIT
   END DO LOOP_I
   IP = IP + 1
END DO
END SUBROUTINE POLY2TRI

END MODULE BOXTETRA_ROUTINES
